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ABSTRACT 


This research considers the theoretical and applied aspects of successive 
approximation techniques for the determination of controls for nonlinear dynamical 
systems. Particular emphasis is placed upon the methods of contraction mappings 
and modified contraction mappings . It is shown that application of the Pontryagin 
principle to the optimal nonlinear regulator problem results in necessary con- 
ditions for optimality in the form of a two point boundary value problem (TPBVP) . 
The TPBVP is represented by an operator equation and functional analytic results 
on the iterative solution of operator equations are applied. The general con- 
vergence theorems are translated and applied to those operators arising from 
the optimal regulation of nonlinear systems. It is shown that simply structured 
matrices and similarity transformations may be used to facilitate the calculation 
of the matrix Green's functions and the evaluation of the convergence criteria. 

A controllability theory based on the integral representation of TPBVP 's, the 
implicit function theorem, and contraction mappings is developed for nonlinear 
dynamical systems. Contraction mappings is theoretically and practically applied 
to a nonlinear control problem with bounded input control, and the Lipschitz 
norm is used to prove convergence for the nondifferentiable operator. A dynamic 
model representing community drug usage is developed and the contraction mappings 
method is used to study the optimal regulation of the nonlinear system. 
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CHAPTER I 


INTRODUCTION 


1.1. Background 

Optimal control theory has experienced an increasing growth of interest in 
the past two decades. Initially motivated by the aerospace effort, optimal control 
theory is now involved in many aspects of general systems engineering. Applica- 
tions range from chemical process control to attempts at managerial and economic 
planning. 

One of the most important and most widely treated problems to date in 
optimal control theory is the so-called "Linear Regulator Problem". Historically, 
this problem arose in Wiener's work concerning stationary time series and linear 
filtering and prediction [Wl] . Under the name "Minimum Integral Squared Error", 
development of this problem was continued through the 1950's by Newton [Nl], 

Booten [B3],and Zadeh [Zl] . Finally utilizing the techniques of modern control 
theory, Kalman [Kl] presented important new aspects of the problem. 

The prominence of this problem is due to two primary factors. First, the 
problem provides a strong link between the classical methods of analytic feedback 
system design via frequency domain methods and the more recent variational approach 
favoring analysis in the time domain [K2], [W2] . Secondly, the problem allows the 
determination of optimal controls in closed form with mathematical ease. (For 
general development and presentation of the problem, see Athans and Falb’ [Al] and 
Lee and Markus [LI]). Finally, a pragmatic motivation for considering the 
problem is the ease with which the quadratic cost criteria can be interpreted 


/ 


1 


physically. Consequently, optimal linear regulation has been extensively 
applied to various systems. For example, the theory has found widespread 
applications in the area of automatic flight control systems. Much of this 
work is based on the significant efforts of Rynaski [R4] , [R5] . Other examples 
of optimal linear regulation are contained in Dyer and McReynolds [D2] . 

However, few systems can adequately be described by a linear dynamic model. 
In particular, increasing effort is now being devoted to the development of 
models representing systems as varied and as complex as urban areas, natural 
resource depletion, management of R and D efforts, and drug usage within a 
community. These models are primarily due to the efforts of Forrester [F3, 
F4,F5,F6] and Roberts [R2] . Along with many engineering systems, these systems 
contain inherent nonlinearities which must be included in any meaningful study. 

In contrast to linear systems, the regulation of nonlinear dynamical systems 
has received limited attention, most of a specialized nature. The primary reason 
for this seems to lie in the fact that nonlinear optimal control problems can 
rarely be solved analytically or, more specifically, in feedback form as for 
linear regulators. As a result, one must often resort to iterative numerical 
techniques for the determination of the optimizing solutions. Consequently, much 
of the analysis regarding regulation of nonlinear systems concerns techniques for 
determining suboptimal feedback controllers. (See for example [Dl] , [G2] , [L3] , 
[PI], [S2] , [Jl], [F8] , and [B5]). Most of these approaches involve the modeling 
of. the nonlinear system as a linear system in some manner. A somewhat different 
approach, not suboptimal, is taken by Brunovsky [B4] and Lukes [L4] . Both of 
these treatments are closely related to the basic hypothesis that the system be 
stabilizable [LI] . Under the assumption of complete controllability, Brunovsky 
approached the problem via Lyapunov functions. Lukes requires the system be 
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stabilizable and then uses Lyapunov-like theory to obtain results for feedback 
controllers . 

The direction of these various approaches is primarily generated by the 
desire for a feedback controller. However, there is a second, more esoteric 
reason, and that is the desire for general results. Unfortunately, the undis- 
cerning application of an algorithm often limits insight into the underlying 
structure of the problem being considered. This loss of general information is 
often due to the fact that practical convergence criteria are few for most of 
the iterative methods used in the solution of optimal control problems. Theoreti- 
cal aspects of these criteria have been investigated by numerous applied mathe- 
maticians (see Kantorovich [K4] and Collatz [C2]). The Russian Kantorovich [K4] 
was one of the first to develop and unify the mathematical theory of iterative 
methods. Using the power of functional analysis methods, he presented conver- 
gence results for such basic iterative schemes as contraction mappings and 
Newton's method. These basic results have been considerably broadened, modernized, 
and made practical by the efforts of Falb and de Jong [FI]. In their book, they 
present the derivation of general convergence criteria for the application of 
various successive approximation methods to the solution of optimal control 
problems . 

1.2. Description of the Problem 

The primary goal of this research is the consideration of the theoretical 
and applied aspects of successive approximation techniques for the solution of 
optimal nonlinear regulator problems. Application of the Pontryagin principle 
to the posed optimization problem results in necessary conditions for optimality 
in the form of a two point boundary value problem (TPBVP), Hence, the central 
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theme of this study shall be the application of successive approximation methods 
to the solution of nonlinear TPBVP's which arise from optimal nonlinear regulation. 
The basic approach to be used is to represent the TPBVP by an operator equation 
and then apply functional analytic results in the iterative solution of operator 
equations . 

In particular, we shall investigate the contraction mappings method and the 
modified contraction mappings method. We have as our first objective the trans- 
lation and application of the general convergence theorems to those operators 
originating in the optimal regulation of a nonlinear system. A second objective 
is the development of techniques to facilitate the evaluation of the convergence 
criteria. Finally, example problems will be solved to demonstrate the usefulness 
of the theory. 

1.3. Synopsis 

A brief summary of the dissertation is as follows: In Chapter 2, the optimal 

regulation of dynamical systems is introduced. In particular, we discuss the 
reduction of optimization problems to two point boundary value problems by means 
of Pontryagin's principle. Results are derived for optimal regulation of linear 
dynamical systems (Section 2.2) and several classes of nonlinear systems (Sec- 
tion 2.3). Optimal system regulation is considered for both unconstrained and 
bounded controls. In Chapter 3, methods of solving two point boundary value 
problems are presented. In particular, the integral equation representation 
of two point boundary value problems is introduced (Section 3.2). The book by 
Falb and de Jong [FI] was used as the main reference for this chapter. The 
integral representation makes it possible to consider the solution of a two 
point boundary value problem as the solution of a corresponding operator equation. 
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A review of Lipschitz norms and derivative norms for the integral operator is 
presented (Section 3.3) and the methods of contraction mappings (Section 3.4) 
and modified contraction mappings (Section 3.5) are introduced. Convergence 
theorems for both methods are presented. Chapter 3 concludes with the application 
of contraction mappings to the solution of two point boundary value problems 
arising in Chapter 2 and the derivation of translated convergence theorems. 

Chapter 4 is devoted to a rather detailed investigation into the calculation 
of the theoretical convergence criteria. Upper bounds are presented for the 
Lipschitz norm and derivative norm (Section 4.2) and various techniques for 
evaluating these bounds are introduced. In particular, the use of simply 
structured matrices (Sections 4.4, 4.5) and similarity transformations (Section 
4.6) are considered. The use of partitioned matrices in these developments 
provides considerable insight into the generic structure of the Green's matrices 
contained within the integral representation. In Chapter 5 the issue of con- 
trollability for nonlinear systems is considered. Specifically, it is shown 
that controllability for linear systems (Section 5.2) and nonlinear systems 
(Section 5.3) may be studied via the integral representation and contraction 
mappings. In Chapter 6 we present numerical examples to illustrate the theoreti- 
cal and practical application of contraction mappings to the regulation and 
control of nonlinear systems. In Chapter 7, a dynamic model is developed for 
a socio-economic system and contraction mappings is used to investigate the 
optimal regulation of this nonlinear system. Finally, in Chapter 8, we summarize 
our results and indicate directions in which future research may be done. We 
conclude with an appendix which gives the actual computer. program (written in 
the FORTRAN language) which was used in the application of contraction mappings 
to the problem discussed in Chapter 7. 
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CHAPTER 2 


OPTIMAL REGULATION OF DYNAMICAL SYSTEMS 


2.1. Introduction 

An optimal control problem is a composite concept consisting of four basic 
elements: (1) a dynamical system, (2) a set of initial states and a set of final 

states, (3) a set of admissible controls, and (4) a cost functional to be minimized. 
The problem consists of finding the admissible control which transfers the state 
of the dynamical system from the set of initial states to the set of final states 
and, in so doing, minimizes the cost functional. In this chapter we discuss the 
optimal regulation of nonlinear systems and the reduction of the optimization 
problem to a TPBVP by means of Pontryagin's principle. 

2.2. Optimal Linear Regulator 

As a preface to the nonlinear system analysis, we shall present the basic 
results for the optimal linear regulator. (For a very thorough treatment of this 
problem see Kleinman [K4]). 

Definition 2.2.1. Linear Dynamical System 

A linear dynamical system is characterized by the following elements: 

(1) A state vector x of dimension n 

(2) A control input vector u of dimension r 

(3) A linear differential equation which describes the evolution of the 
system in time, i.e., 

i(t) = A(t) x(t) + B (t) u(t) 2.2.2 

where A(t) is an n*n matrix and B(t) is an n*r matrix. 
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Now given an initial state, x(t 0 )= x^, and assuming the control u[t) is not 
constrained, the optimal linear regulator problem is then to determine the control 
u(t) which minimizes the quadratic cost function 

J(u) = i<x(T), Kx(T) > + \ J [<x(t), Q(t) x(t)> + <u(t), R(t) u(t)>]dt 

0 2.2.3 

where 

The terminal time T is specified 2.2.4 

K ds a constant nxn positive semidefinite matrix 

, Q(t) is an nxn positive semidefinite matrix 

✓ 

R(t) is an r*r positive definite matrix 
aqd k and Q(t) are not both identically zero. 

Application of the- minimum principle to the optimization problem posed above yields 
necessary conditions for optimality in the form of the 2n x 2n canonical system 
o.f equations 

■ *i r 

A(t) 

|-QCt) 


x(ty 


p(t) 

■ 


)R _1 (t) 1 
-A'(t) 



x(t) 

. 

p(t) 


2.2.5 


subject to the boundary conditions 

t 

X fV' = X 0 

■p.m = K x(T) . 2.2.6 

O 

The H-minimal control for'tE[tg, T] is then given by 

u(t) = -R _i Ct) B ' (t) p(t). 2.2.7 

The boundary conditions specified by eq(2.2.6) may be expressed more compactly as 


0-3 rx(t 0 )-l + ro o-|rx(T)l rx 0 - 

-6 Oj Lp(t 0 ) J + L-K I J L P CT)J Lo. 


2 . 2.8 


where I is the nxji, identity matrix and 0 is the nxn zero matrix. This form of 

- ' - 

expressing boundary conditions will become important in the sequel. The TPBVP 
arising from the linear optimal regulator problem may then be put into the form 
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y(t) = S(t) y (t) 

My(tg) + Ny(T) = c 
where y is the 2n composite vector 
'x(t)' 

y(t) = 

p(t) 

- 

S(t) is the 2n x 2n matrix 
T Aft) 


,-l. 


set) = 


-Q(t) 


-B(t]R (t) B ' (t) 
-A'(t) 


N are the 

2n 

x 2n boundary 

value 

r 1 

°-| 

r 0 

°-| 

M = 


N = 


L 0 

0 J 

L-k 

iJ 


and c is the 2n constant matrix 




2.2.9 


2 . 2.10 


2.2.11 


2. 2. 1'2 


"2.2.13 


In many physical situations, the input control u(t) may not take on all values 
As an introduction to systems with bounded control, let us now suppose the input 
control to the linear dynamical system is constrained in magnitude by the relation 
| Uj CO| II j = l,...,r . ' • • 2.2.14 

Then given an initial state for the linear dynamical system, the optimal linear 
regulator problem is to determine an admissible control u(t)6fi which minimizes 
the quadratic cost functional given in (2.2.3). 

If is shown in [Al] that the necessary conditions for optimality reduce to 
the 2n x 2n canonical system of equations 

x (t) = A(t) x (t ) - B(t) SAT (R _1 (t) B'(t) p(t)} . 2.2.15 

p(t) = Q(t) x(t) - A' (t) p(t) ’ • 

subject to the boundary conditions 

x (t Q ) = x 0 2.2.16 

P(T) = K x(T) , 
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where the SAT function is defined as 


i , y > i 

SAT{y } = { y , |y| <_ 1 2.2.17 

•1 , y < -1 

It is seen this system of 2n differential equations is not linear. The necessary 
conditions thus reduce to a nonlinear TPBVP of the form 

y(t) = set) y(t) + f(y(t)) 2.2.18 

My(t Q ) + Ny(T) = c 
where y is the composite 2n vector 

r *(t)n 


y(t) = 

L p(t). 

S(t) is the 2n x 2n matrix 

A(t) 0 
-Q(t) -A'(t). 

f Cy ) ) is the 2n vector function 


S(t) 


rA(t) 0 -| 

L-Q(t) -A'(t)J 


f(y(t)) 


■[ 


-B (t) SAT{R _1 (t) B' (t)p(t) } 


] 


M and N are 2n x 2n matrices and c is the 2n vector 


r 1 °i » = r 0 "i .. r»i 

- 0 0 J L -K I J Lo . 


0 0 


2.2.19 


2 . 2.20 


2 . 2.21 


2 . 2.22 


This example illustrates a nonlinear TPBVP arising from the optimization of a 
linear system. We shall now consider the optimization of nonlinear systems 
and the forms of the resultant TPBVP' s. 
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2.3. Optimal Regulation of Nonlinear Systems 

In this section we shall consider the control of several classes of nonlinear 
systems subject to the quadratic cost functional given in (2.2.3). Our aim in 
this section is to reduce the necessary conditions for optimality to two point 
boundary value problems. 

Example 2.3.1. 

Many nonlinear systems contain nonlinearities involving only the state 
variables. Hence, rather than initially considering the most general formulation, 
we shall first consider the class described by the differential equation 


x(t) = A(t) x (t ) + B (t) u(t) + ij,(x(t)) 


2.3.2 


where we assume iji(x(*)) and (3 i(i/3x) (x(- ) ) are continuous on R n . We shall initially 
consider the control to be unconstrained, i.e., u€Q = R n . Again we shall 
consider the quadratic cost functional 


J(u) = ±-<x(T),Kx(T)> + 


1 

2 



<x(t), Q(t) x(t)> + <u(t), R(t) u 


(t)>]dt 

2.3.3 


subject to the assumptions of (2.2.4). Application of the minimum principle 
yields necessary conditions for optimality in the form of the 2n * 2n canonical 
system of equations 


X(t) 


'A(t) 

~B (t) R 1 (t) B ' (t) 

rx(tr 


*(x(t)) 

p(t)_ 


Q(t) 

-A’(t) ■ 

LpW. 


- (3*/3x)'(x(t)) p(t). 


subject to the boundary conditions 


I 0 


X( V 

+ 

0 

O' 


• X CO- 

= 

x o 

0 0 


pcto). 


-K 

I. 


PCO. 


0 


The H-minimal control for t£[tg,T] is then given by 

u (t ) = -R _1 (t) B'(t) p(t). 2.3.6 


It is often advantageous to standardize the time interval over which the TPBVP is 
defined. This standardization is accomplished by the introduction of a new 
variable. Let (see Long [L2]) 
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Here s is the new variable which varies between 0 and 1. In most cases we may 
take tg = b = 0. In terms of s and a, the TPBVP then becomes 


' X(s) ■ 

= a 

A (as) 

-B(as)R ^(as)B'(as) 

x(s)‘ 

+ a 

Kx(s)) 

. P(s) . 


. -Q(as) 

-A' (as) 

p(s). 


- (3ij»/3x) * ( x(s ) ) p (s ) 


2.3.8 


with 


I 0 

x (0) 


0 

0 

'x(l) 


x o 

0 0 

P(0) 


-K 

I 

,P(1). 


0 . 


where (') indicates differentiation with respect to s. In the sequel, the 
TPBVP's which shall be considered will generally be normalized in this fashion. 
Example 2.3.10. 

As an illustration of the ideas presented in Example (2.3.1), let us consider 
the driven, second order nonlinear oscillator studied by Van der Pol. We have 
the system given as 


x x (t) = x 2 (t) 2.3.11 

* 2 (t) = + «Cl-*J(t))x 2 (t) + u(t). 


or in vector-matrix form as 


X 1 


0 1' 


x l' 


■ 0 ' 


0 

1 

= 



+ 


u + 

2 

X 2 


-1 0 . 


X 2‘ 


. 1 . 


e(l-x 1 )x 2 . 


The optimization problem to be considered is that of minimizing the cost 
functional 

T 

J = j f (xf(t) * x 2 (t) + u 2 (t))dt 
0 


2.3.13 



subject to the boundary conditions 


x^O) = x Q , x^T) unspecified 
x 2 (0) = 0 , x 2 (T) unspecified. 

From eq (2.3.4), we have the 2n x 2n canonical system 


**1 


0 

1 

0 

0 

X 1 


0 

*2 


-1 

0 

0 

1 

X 2 

+ 

e(l-x^)x 2 

Pi 


-1 

0 

0 

1 

Pi 


2ex 1 x 2 P2 

P 2 


0 

-1 

-1 

0 

P 2 


-e(l-x^)p 2 






J 





subject to the boundary conditions 


*1 

0 

0 

o' 


X x (0) 


0 

0 

0 

0 


XjU) 


X 0 

0 

1 

0 

0 


x 2 (0) 


0 

0 

0 

0 


x 2 (i) 


0 

0 

0 

0 

0 


Pj(0) 


0 

0 

1 

0 


p^l) 


0 

0 

0 

0 

0 


! > 

o 

V > 

CM 

a. 


p 

0 

0 

1 


,p 2 a). 


0 . 


or, in the more compact form, 

y(t) = Sy(t') + f (y(t)) 
My(0) + Ny(l) = c. 


2.3.14 


2.3.15 


2.3.17 


2.3.18 


In the sequel, this example will reappear as we consider the iterative solution 
of TPBVP's of the form (2.3.18). 

The class of systems studied in Example 2.3.1 will now be reconsidered with 
a magnitude constraint upon the control. 

Example 2.3.19. 

Let us now consider the regulation of the previous system 


*(t) = A(t) x(t) + B (t) u(t) + ij>(x(t)) 


2.3.20 
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where the input control vector is constrained in magnitude by 


[iij (•) | < 1 , j =1 , . . . ,r . 2.3.21 

The cost functional is again given by (2.3.3). Application of the minimum 
principle yields the 2n x 2n system of canonical equations as 


• X(t)1 


■ A(t) 0 

x(t)' 


'♦(x(t)) - B (t)SAT{ R _1 (t) B ' (t)p (t) } 

. p(t) . 


-Q(t) -A' (t) 

P(t). 


- (a*/a^x(t)) P (t) 


2.3.22 

subject to the boundary conditions 


i 

o' 

/ — * 
O 

v — / 
X 

+ 

' 0 

0 ' 

x(T) 

_ 

x o 

0 

0 

l pCTq) J 


_ -K 

I . 

. P(T). 


0 


where the SAT function is specified in (2.2.16). For this system, the H-minimal 
control is given as 

u (t) = -SAT{R _1 (t) B'(t) p(t)} , te[t 0 ,T], 2.3.24 

Example 2.3.25 

In the previous example, we discussed the laTge class of nonlinear systems in 
which the nonlinearity is a function of only the state variable. Let us now 
consider the more general system described by the differential equation 


x = A(t)x + B(t)u + i^(x,u) 


2.3.26 


where A(t) is an n x n matrix, B(t) is an n x r matrix, u is an unconstrained 
r-vector, and iJj(x,u) and (3i|//8x) (x(- ) ,u(- ) ) are continuous in R n x R n . The 
system is subject to the quadratic cost criteria given as 


T 

J = j <x(T) ,Kx(T)> + j f [<x(t), Q(t)x(t)> + <u(t), R(t)u(t) > ]dt 

2.3.27 
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under the assumptions of (2.2.4). 

Following the Pontryagin minimum principle, the Hamiltonian for the. 
optimization problem posed above is given as 

H = i <x (t) , Q(t) x(t)> + y<u(t), R(t)u(t)> + <A(t) x(t), p(t)> 

2.3.28 

+ <B(t) uCt), p(t)>+ <iHx(t) , u(t)), p(t)> . 

Formally applying the Pontryagin principle, the costate vector is then described 
by the differential equation 

p(t) = - Q(t) x(t) - A'(t) p(t) - (3^/3x)'(x(t),u(t))p(t). 2.3.29 

Along the optimal trajectory we must have u*(t) minimizing the Hamiltonian, i. e. , 

H(x* (t) , p*(t), u*(t),t) < H(x*(t),p*(t),u,t) 2.3.30 

for all admissible w and where (•)* denotes optimal trajectories. If the 
Hamiltonian is normal [Al], the minimization equation (2.3.30) may be solved 
for the H-minimal u in terms of x,p, and t, i.e., 

u = ?(x,p,t) . 2.3.31 

Now using (2.3.31) we define 

l Hx,?(x,p)) = <Hx,p) 2.3.32 ' - 

and 

(3i|i/3x) (x,5(x,p)) = Z(x,p) 2.3.33 

where <f> is an n vector function and Z is an n x n matrix function. 

xl A(t) 0 I [x I [B(t)£(x,p) + Kx,p) 

= + 
pJ l -Q(t) -A’ (t) J [p J - Z'(x,p)p 

subject to the boundary conditions 
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2.3.35 


' I 

o' 

/ N 

O 

X 

i . ... 

+ 

0 

o' 

x (T) 

= 

V 

.0 

0 . 

1 PCt 0 ). 


-K 

I 

P(T). 


o . 


These results may be applied to various forms of ijj(x,p). In the following 
example we consider one such form. 

Example 2.3.36 

Consider the system described by 

x = A(t)x + B(t)u + D(x)u 2.3.36 


where A(t) is an n x n matrix, B(t) is an n x r matrix, D(x) is an n x r matrix, 
D. . (x) and (3D. ./3x, ) (x(-)) are continuous in R n , and u(-) is an unconstrained 

1J 1J K 

r vectof. Consider system (2.3.36) subject to the cost functional (2.3.27) and 
the initial condition x(t Q ) = x 0 - 

Define the vector function 5(x,p,t) to contain the elements 


5^Cx,p,t) = <^p, - (3D/3x^) (x) R' 1 (t)[B'(t) + D'(x)]^> 
and define the matrix c(x,t) as 


2.3.37 


C(x,t) =-B(t) R 1 (t) D'Ct) - D(x) R' 1 (t)[B , (t) + D'(x)J. • 2.3.38 

Using the results of Example 2.3.25 and (2.3.37), (2-3.38), the 2n x 2n canonical 
system of equations is given as 


it 

_ 

A(t) 

-B(t)R'l(t)B’ (t)' 

x‘ 

+ 

C (x,t)p ■ 

P 


-QCt) 

-A(t) 

P. 


5(x,p,t) 

I 

O' 

xC V 

.1° °1[ X(T, 1 

* X 

0 


0 

0 

IpCV- 

1 

►— « 

' — > 
H 
' — ' 

. 0 




The various two point boundary value problems presented in the previous 
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examples can very rarely be solved analytically. Thus we shall investigate 
successive approximation techniques for the solution of TPBVP's in the next 
chapter. 
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CHAPTER 3 


METHODS OF SOLVING TPBVP's 


3.1. Introduction 

In the analysis of optimal control problems, the necessary conditions for 
optimality are often in a form which may be reduced to a TPBVP of the form 

y(t) = F(y.t), gfy(O)) + h(y(l)) = c. 3.1.1 

In particular, we presented in Chapter 2 various TPBVP's which originate in the 
optimal regulation of certain classes of nonlinear systems. We shall now illustrate 
that under certain conditions, such TPBVP's may be represented by operator equations 
of the form 

y=T(y). 3.1.2 

Then, following the lead of Falb and deJong [FI], we shall investigate the applica- 
tion of successive approximation techniques to the iterative solution of these 
operator equations . 

3.2 Representation of TPBVP's 

In this section we consider the (normalized) two point boundary value problem 
y(t) = F(y,t) 7 g(y (0)) + h (y (1) ) = c 3.2.1 

p 

where G, g, and h are vector valued functions and c is an element of R . We 
shall first review some results relating to the development of equivalent integral 
equation representations of the TPBVP (3. 2 . 1) . Most results in this section. 
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come from Falb and de Jong [FI]. Since linear TPBVP's will play an important 
role in the integral equation representations, we begin our discussion with a 
consideration of linear TPBVP's. 

Consider the linear TPBVP 

y(t) = V (t)y (t) + f(t) , My (0) + Ny(l) = c 3.2.2 

where V(t), M, and N are p k p matrices, and f(t) and d axe p vectors. We present 
the following theorem on the existence of a solution of equation (3.2.2). 

Theorem 3.2.3 

Suppose that the functions V(t) and f(t) satisfy appropriate smoothness and 

V V 

boundedness conditions and det[M + N$ (1,0)] ^ 0 where $ (t,s) is the fundamental 
matrix of y = Vy. Then (3.2.2) has a unique solution y(t) on [0,1] which can 
be written in the form 

1 

y (t) = H (t) c + J G J (t,s)f(s)ds 3.2.4 

0 

where the Green's matrices H and G are given by 

H(t) = $ V (t ,0) [M+N$ V (1 ,0) ] -1 3.2.5 

and 

f «J> V (t,0) [M+N# V (l,0)] -1 M* V (0,s) , 0< s< t 3.2.6 

V V -1 V 

-# (t,0)[M+N4> (1,0)] N<J (1 ,s) , t < s^l 

for all t,s in [0,1] . 

Proof: (See [FI] for proof of theorem and technical conditions specified for 

V(t) and f (t)) . 
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The requirement in Theorem 3.2.3 that det [M+N^ (1 ,0) ] ^ 0 is crucial to the 
integral representation of TPBVP's. We therefore make the following definition. 
Definition 3.2.7 

Let V, M, N be p x p matrices. Then J = {V(t),M,N} is called a boundary 

compatible set if and only if V(t) satisfies certain technical conditions and 
V V 

det[M+N$ (1,0)] i 0 where 4> (t,s) is the fundamental matrix solution of 
y = V(t) y. 

In the sequel we shall often be given two boundary related matrices M and N and 
will be required to determine a matrix V(t) so that the set J = {V(t),M,N} is 
boundary compatible. In the next lemma we give necessary and sufficient 
conditions for the existence of a matrix V(t) which is boundary compatible with 
the prescribed matrices M and N. 

Lemma 3.2.8 

Let M and N be p * p matrices. A necessary and sufficient condition that 
there be a V(t) with J = {V(t),M,N} boundary compatible is that the p x 2p matrix 
[M N] have full rank p. 

Proof: (See [FI].) 

Theorem 3.2.3 and Lemma 3.2.8 form the basis for the integral equation representa- 
tion of nonlinear TPBVP's of the form (3.2.1). We now have the following. 

Theorem 3.2.9 

Suppose that F(y,t) satisfies certain technical conditions and J = (V(t),M,N) 
is a boundary compatible set of dimension p. Then the boundary value problem 

Y = F(y,t) 7 g(y (0)) + h (y (1) ) = c 3.2.10 

has the equivalent representation 
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3.2.11 


y(t) = H (t){c - g(y(0) - h(y(l)) + My(0) + Ny [1) } 
1 


f Ct ,s) {F(y (s) ,s) - V (s)y(s) }ds 


where the Green's functions H^(t) and G^(t,s) are given by 


H J (t) = 4> V (t ,0) [M + N$ V {1 ,0) ] " 1 


and 


G J (t,s) 


<i> V (t,s)[M + N$> V (l,0)] _1 M$ V (0,s) , 0 < s < t 


-<f V Ct,s)[M + N* V (l,0)] _1 N* V (l,s) , t < s < 1 


3.2.12 


3.2.13 


V 

where $ (t,s) is the fundamental matrix of the linear system y = V(t)y. 

Proof: (See [FI] for complete conditions assumed for F(y,t) and a proof of 
the theorem.) 

Theorem 3.2.9 presents an integral equation representation for TPBVP's of 
the form (3.2.1). It is now a simple matter to demonstrate that solving (3.2.1) 
is equivalent to solving a certain fixed point problem in an appropriate Banach 
space. In particular, assuming that the conditions of the previous theorem are 
satisfied, we can define a mapping of the Banach space Y = ^([0,1],R P ) into 
itself by setting 


T J (y) = H J (t){p-g(y(0)) - h(y(l)) + My(0) + Ny(l)} 

1 

+ J G J (t,s)(F(y(s) ,s) - V (s)y (s) }ds . 

0 


3.2.14 


Then, (3.2.11) is equivalent to the fixed point problem 


/ = T J (y) 


3.2.15 
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on ^ ( [0 , 1 ] , rP) . The operator equation (3.2.14) can now be solved by successive 
approximation iterative techniques as presented by Kantorovich [K4] and particu- 
larly Falb and de Jong [FI]. 


3.3. Frechet Derivatives and Lipschitz Norms 

In the discussion of successive approximation iterative techniques, we shall 
require an expression for the Frechet derivative or Lipschitz norm of the operator 
T^. In this section we shall present a brief treatment of these concepts. (Again, 
many of these basic results are from Falb [FI].) Let us begin with the following 
definition . 

Definition 3.3.1. . . 

Let Y be a Banach space with |) • || as norm. Let f! be a closed subset of Y 
and let T map Y into Y. The Lipschitz norm of T on ft, in symbols: 0T 0 ^ , is 
given by 

S T lfi = u*v*ft t II T (u) - T (v) || / |l u-v II } . 3.3.2 


If T is Frechet differentiable on ft, then derivative norm of T on ft, in symbols: 
0 T fl ^ is given by 

I 'V" • 3 - 3 ' 3 

We shall now compute expressions for (T 1 and (T ^)" . We have 

(T y J )’(u) = H J (t)([M-(3g/3y)(y(0))]u(0) ♦ [N- (3h/3y) (y (1) ) ]u(l) } 

1 

+. f G J (t,s)((3F/ay)(y(s)) - V(s))u(s)ds 3.3.4 
0 


and 
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/ V 

(T^j) "(u, v) = H J (t)p[(3/3y i )C-3g/8y)](y(0))u.(0)v(0) 


i=l 


2 tC3/3y i ) C-3h/3y)Ky(l))u. (l)v(l) | 


i=l 


* / ° J(t - 


S) Y1 [(3/3y i )(3F/3y)](y(s))u.(s)vCs) ds 3.3.5 
' i=l 


provided the indicated partial derivatives exist. When evaluating convergence 


criteria, we shall require estimates, say for example of the norm of the operator 
(T ' . There are of course several expressions for calculating or estimating 


(T y J ) 'll 


Since the more accurate expressions are difficult to evaluate in 
practice, we shall present a coarse estimate that is more amenable to future 
applications. We recall first of all that if v(>)£ ^([0,1], R^) , then 


II v ( - ) B =■ sup sup |v. (t)| 3.3.6 

icP tt[0,l] 1 

is the norm of v(-) where P = {l,...,p} and v i (‘) is the ith component of v(-)* 

Noting that ll(T J ) 1 1| = sup j || (T J )'ull! and letting H^ft) = [H^.(t)], 
y Hull <1 1 u 1 1J 

G J (t,s) = [G^j (t, s )] , M = [m^ k ], N = [n , V(s) = [v^ k (s)], we have as a coarse 

estimate 


|J(T J )'|! = sup |||(T J )’u||| 

3 m ii < y > 

! P P 

£ (l H ij( t )l) / t k - (ag J -/ay k )(y(o))| 
j=l ' k=l 

+ l n jk ■ ( 8 h j / 8 y k )(yCD)l>) 3.3.7 
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+ ( f l G ij t t » s ^ ds )'( su p | ^ I ( 8 F j/ 3 y k ) Cy C s ) , s ) ' v jk (s) l{)| 

j=l 0 s k=l 

Expression (3. 3. 7). will become quite important in the sequel. One of our primary 
objectives shall be determining techniques for easily estimating this expression. 

In some cases, the smoothness conditions required to obtain Frechet deriva- 
tives are too strong. As an example, we have the nonlinearity containing the 
SAT function in equation (2.2.15). This fact does not imply that successive 
approximating techniques may not be applied to the iterative solution of the 
operator equation. It simply means we have lost one method of evaluating 
convergence criteria. Hence, under somewhat weaker smoothness conditions, we 
shall compute the Lipschitz norm of the operator T J (y). 

We have the following result from Falb [FI] . 

Lemma 3.3.8 

Let S be a bounded open set in £([0,1], R* 5 ) and let D be an open set in 

rP containing the range of S. Suppose that (i) K(t,y,s) is a map of 

[0,1] x D x [0,1] into D which satisfies certain technical conditions, and (ii) 

there is an integrable function m(t,s) of s with sup / m(t,s)ds = u < “such 

t •'O 

that I K(t,y,s) II <_ m(t,s) and ||K(t,y 1 ,s) - K(t,y 2 ,s)|| <_ m(t,s) ||y 1 -y 2 || on 

[0,1] x D x [0,1], Then the mapping T given by 

1 

T(u)(t) = / K(t,u(s) ,s)ds 

0 

maps £([o,i], R^) into £([0,1], R^) and the Lipschitz norm, 3 T 0 g, satisfies 

I T 0 <_ y. 3.3.9 

Proof: (See [FI] for proof of theorem and specific conditions on K.) 
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Corrollary 3.3.10 

Suppose that the function 

K(t,y,s) = G J (t,s){F(y,s) - V{s)y} 
satisfies the conditions of Lemma 3.3.8 and that 

ii g(y 2 ) - g(y 2 ) ii I'Wjii y 2 -y 2 it 311(1 11 h (y 1 )- h (y 2 ^> 11 1 v 2 ii yj-y^l 3 - 3 - 11 

Let 

a = max {y, || H J (0 II V l , » ^ ( • ) || v 2 , || H J ( • )M II ., J H J ( ■ )N II } . 3.3.12 

Then BT J 0s 4= « 3.3.13 

This result will prove useful in particular when investigating regulators with 

bounded input controls. 

3.4. Contraction -Mappings Method 

Contraction mappings (or Picard's method, [ P2 ] ) is well known in the mathemat 
ical literature and has long been a standard approach for proving existence and 
uniqueness properties for ordinary differential equations. (See for example 
Coddington and Levinson [Cl], specifically Section 1.3 entitled "The Method of 
Successive Approximations.") To formalize our discussion of this technique, let 
us begin with the following definition. 

Definition 3.4.1 

Let Y be a topological space and let T map Y into itself. Let y^ be 
an element of Y. The sequence (y (■)} generated by the algorithm 

y n+ i = T Cy n ) n= 0,1,2,... 3.4.2 

is called a contraction mapping or CM sequence for T based on y^. 

The following theorem is central to our future discussions concerning the 
contraction mappings method. 
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Theorem 3.4.3 


Let Y be a Banach space and let S = S(y Q ,r) be the closed sphere in Y with 
center y^ and radius r. Let T map Y into Y and suppose that (i) T is defined 
on SCy^.r), and (ii) there are real numbers n and a with n 0 and 0 _< a < 1 
such that 


f y x - y 0 I* I n 

Blfig-^cicl or lT§ , g<_a<l 


3.4.4 


3.4.5 

3.4.6 


where y^ = TCy^). Then the CM sequence {y^} for T based on y^ converges 
unique fixed point y* of T in S and the rate of convergence is given by 


UU LUG 


y* - y r 


— l-a 


y n ' VlH 


3.4.7 


Proof: (See [F ]). 

Let us now consider the application of this theorem to operator equations of the 
form 


y(t) = T J (y)(t) = H J (t) {c-g(y (0) ) - h(y(l) * M y(0) + Ny(l)} 

1 3 
+ f C J (t,s){F(y(s) ,s) - V(s)y(s)}ds 
0 

where J = {V(t),M,N} is a boundary compatible set. Following the contraction 
mapping prescription, we select an initial element YqCO in ^([0,1] > R ) and 
successively generate a CM sequence ,{y (•)) f° r based on ( * ) by means of 
the algorithm 

Y , = T J (y ) 3 

■'n+1 w n 
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or equivalently, by 


y n+ i CtD = H J (t){c-gCy n (0)) - h(y n (0)) + My n (0) + Ny (1)} 

3.4.10 

+ f G J (t,s){F(y n (s),s) - V(s)y n (s)}ds. 

Since we know y (•) at each successive step, we can write (3.4.10) in the form 

y n + l (t) = + f G J (t,s)f n (s)ds 3.4.11 

where 

c n = c - g(y n (0)) - h(y n (l)) + My n (0) + Ny n (l) 3.4.12 

and 

f n (s) = F(y n (s)) - V(s) y n (s). 3.4.13 

Hence, it is seen from (3.4.11) and our results on linear TPBVP (eq. 3.2.4) that 
the method of contraction mappings when applied to (3.4.8) essentially amounts 
to the successive solution of the linear TPBVP's (3.4.11). 

If the partial derivatives of (3.3.) exist, we then have the following. 
Theorem 3.4.14. 

Let y Q (-) be an element of £([0,1],R P ). and let S = S(y Q ,r). Suppose that 
(i) J = (V(t),M,n) is a boundary compatible set for which 

y = F(y(t),t) g(y (0)) + h(y(l)) = c 3.4.15 

is differentiable on S, and (ii) there are real numbers n and a with n > 0 and 
0 < a < 1 such that 
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3.4.16 


*T J (y 0 ) - y 0 II = sup sup {|T J (y ).(t) - y„ . (t) | } <_ n 

0 0 i te [ 0 , 1 ] u 0,1 

sup { II (T y J ) ' II } <_ a 3.4.17 

y * s 

n < r - - 3.4.18 

1-a — 

Then the CM sequence {y ( * ) } for the TPBVP based on y Q and J converges uniformly 
to the unique solution y*(-) of (3. 4. IS) in S and the rate of convergence is 
given by 

II y*CO - y n (0 | i || Xi (*) - y 0 (O II- . 3.4.19 

Proof: Simply apply Theorem 3.4.3. 

It should be noticed that if the TPBVP of interest is not differentiable, 
but a Lipschitz norm can be obtained, then (3.4.17) is simply replaced by 

§ < a . 3.4.20 

We shall use (3.4.20) in the investigation of optimal regulators with bounded 
control . 

At this point we shall make a few general comments concerning our representa- 
tion of TPBVP 's and, in partiuclar, the role of the boundary compatible set. 

J = { V,M,N} . From Theorem 3.4.14, we see that the convergence rate factor, a, 
is determined by the Frechet derivative of the operator T^(y). In particular, 
from equation 3.3.6 we have an estimate for this norm given as 
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B(T J )'II = 

y 



4 ± fgP s \?p|X! ( I ^ t)} I } -(E ! |m jk"(3f L ) (y(0)) l+|n j k "(^ 1 ) (y(1))| |) 

( j=l k=l k k 

+ s(i |G ±j (t » s) l ds )* ( s gp jZlf^ 1 ) (y(8),s)-v (s)|j)J 

j=l 0 *k = l k JK ; ; 


3.4.21 


For convergence purposes we wish to make this quantity as small as possible, 
and in this light, we shall discuss the choice of J = (V(t), M,N}. All of the 
TPBVP's obtained in Chapter 2 have linear boundary conditions of the form 
Ky(0) + Ly(l) = c. From this we shall clearly choose M and N to equal the 
linear boundary conditions of the TPBVP, thus eliminating the first terms in 
(3.4.21). We then have the simplified expression 



Consideration of this expression allows us to deduce that if y Q is a good initial 
estimate of the solution, then it is often effective to choose V(s) close to 
(3F/3y) (y Q (s) ,s) . In fact, for V(s) = (3F/9y) (y^ (s) ,s) , the iterative method 
is known as the "modified Newton's method." However, a general choice such as 
this for 
the term 

criteria. In the next section we shall consider a technique which is often 
useful for evaluating convergence criteria. 


ie V matrix usually precludes any attempt at calculating or estimating 


,s)]ds , thus preventing an easy estimation of the convergence 
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3.5. Modified Contration Mappings 

In some situations, the direct application of the contraction mappings 
method does not lead to a convergent sequence of approximations. However, 
it is frequently possible to modify T in such a way as to lead to a convergent 
sequence of approximations. We consider the following. 

Lemma 3.5.1. 

Let T and U be maps of Y into Y. Suppose that I - U is invertible and let 
P. be the map of Y into itself given by 

P (y) = [I-U] _1 [T(y) - U (y) ] . 3.5.2 

Then y*(0 is a fixed point of T if and only if y*(-) is a fixed point of P. 
Proof: (See [FI]). 

We shall then consider the selection of an initial approximate solution y^ and 
the generation of a sequence {y^} by the algorithm 

y n + l = P(y n 3 = [I_U] 1 [T( V ' UCy n )] ‘ 3 - 5 - 3 

We shall call this algorithm the modified contraction mappings method. It 
should be noted that the modified contraction mapping sequence for T based on 
y^ and U coincides with the contraction mapping sequence for P based on y^. 

Hence we may translate the results on contraction mappings into theorems for 
modified contraction mappings. The primary theorem is given as follows. 

Theorem 3.5.4. 

If U is a linear operator with I-U invertible, if T is differentiable on S, 
and if there are real numbers q,a with q > 0 and 0 a < 1 such that 


31 



3.5.5 


II Y x - y Q I in 

sug { II [I-U] 1[T. . - U] II } <_ ci 

yes lyJ 


3.5.6 

3.5.7 


then the modified contraction mappings sequence {y } converges to the unique 
fixed point y* of T and S and the rate of convergence is given by 


y n y n-l A — 1-a 


y l ' y o 


Proof: Apply Theorem 3.4.3. 

The importance of these results lies in the fact that they extend the range of 
applicability of the contraction mapping method to fixed point problems for 
operators T that are not contraction mappings. In other words, the basic 
contraction mapping criteria 

sup { |l (T J ) ' II } < a < 1 3.5.9 

y£S y 

is replaced by the condition that the Frechet derivative satisfies 

sup { l| [I-U]' 1 [T J )'-U] II } < a < 1. 3.5.10 

y £ S y 

A second possibility is to replace the single norm in (3.5.10) by a product 
of two norms so that 


sup { 6 [I-U] _1 || • || [(T J )’-U] II } <_ a < 1. 3.5.11 

y£S y 


This formulation offers the possible advantage of easier evaluation, but also 
results in less sharp convergence conditions. 

We shall now specify the form of linear operator U that will be used in 
the modified contraction mappings algorithm. The following lemma involves 
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the relation between the operators and for different boundary compatible 
sets J = {V(t),M,N} and J = {W(t),K,L}. 

Lemma 3.5.12. 

Let J = {V(t),M,N} and J = {W(t),K,L} be boundary compatible sets. Let 
F(y,t) be continuous in y for each fixed t and measurable in t for each fixed y 
with |F(y,t)|| <_m(t), m(t) integrable. Let T be the linear manifold of 
absolutely continuous functions in {£(.[0,1],$.^). Let be the operator 

given by 

ld(y)(t) = H J (t) (-Ky(O) - Ly (1) + My(0) + NyCl)} 

3.5.13 

+ J G (t ,s)(W(s)y(s) - V (s)y (s) } ds 
0 

for yCO in £([0,1], R P ). Then (i) maps g ([0,1] ,R P ) into g ([0,1] ,R P ) and 
r into r ; (ii) the operator I-U^ L has a bounded linear inverse on r with 

[i-u^r 1 y = [i-v^]y 3.5.14 


for y in T and 

y = H J (t){-My(0) - Ny(l) + Ky(0) + Ly(l)} 

♦/ G J (t»s){V(s)y(s) - W(s)y(s)}ds 

(iii) if yCO is in ^([0,1], R P ), then 

T'V) = [i~u^l ] _1 [T J (y) - u^ L (y)] 


3.5.15 


3.5.16 


and (iv) under the differentiability assumptions 

<V’ ) ' - [I-^tT 1 [ O'/)' - < L J- 3.5.17 

Proof: (See [FI]). 
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We shall limit our future discussions to operators U = of the form given 
by (3.5.13). It then follows from Lemma 3.5.12 that the modified contraction 
mapping method when applied to the equation y = T^(y) with modifying operator 
U = lT,L , is equivalent to the contraction mapping method applied to the equation 

K.L 

(V 

y = T J (y). 

The importantce of this point will become clearer as we develop techniques 
for estimating (T ^) ' . We shall now indicate the approach that will be considered. 
Suppose that J is a boundary compatible set for which the corresponding Green's 
matrices are easy to evaluate and estimate. Then if ||u^ j) <_ q < 1 so that 
|[ U KL^ ^ II — Fq » we can obtain 311 estimate of (3.5.17) which involves only 
the Green's matrices corresponding to J. This advantage may well offset the loss 
of accuracy resulting from using (3.5.11). We now have the following. 

Theorem 3.5.18 

Let y^f-) be an element of ^([0,1], R^) and let S = Sfy^.r). Suppose that 
(i) J = (V (t) ,M,N) is a boundary compatible set for which 

y = F(y,t) g(y (0)) + h(y(l)) = c 3.5.19 

is differentiable on S; (ii) J = {U(t),K,L} is a boundary compatible set; and 
(iii) there are real numbers n,q,8, and a with n 0, 0 <_ q < 1, g >_ 0, and 
a = g/(l-q) < 1 such that 

|T J (y Q ) - y Q If = supsup (|T J (y 0 ) i (t) - y Q i (t)|} <_ n 3.5.20 

it * 

llUj^lUq 3.5.21 

sup { II (T J ) ' - IfL il } < 3 3.5.22 

yes y KL - 

— n < r. 3.5.23 

1-a — 
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Then the MCM sequence (y n (*)} for based on XqC 1 ) and converges uniformly 
to-the unique solution y*(0 of (3.5.19) in S and the rate of convergence is 
given by 

n 

II y*(0 - y n (0 il i y— 1 y 1 (0 - y 0 (O II . 3.5.24 

Proof: Apply Theorem 3.5.4. 

In order to illuminate the preceding discussion, let us consider an example 
utilizing the' previous concepts. 

Example 3.5.25. 

Let us consider the iterative solution of the differentiable TPBVP given as 

y(t) = F(y,t) Ky(0) + Ly(l) = c. 3.5.26 

We shall discuss the choice of the boundary compatible set J = {W(t),M,N} to be 
used in the integral representation of the TPBVP. Since the boundary conditions 
of (3.5.26) are linear, we shall choose M = K and N = L. Let us suppose that 
y 0 (t) is a good initial estimate for the solution of (3.5.26). Then as indicated, 
let us choose W(t) as 

W(t) = (3F/3y)(y 0 (t)), 3.5.27 

assuming this choice of J = (W(t),K,L) is boundary compatible. However, this 

general time varying choice for W(t) makes it extremely difficult, if not 

W 

impossible, to analytically calculate the fundamental matrix 4 (t,s) and the 
Green's functions. 

Let us now decompose the W(t) matrix as 

W(t) = V + 6V(t) 3.5.28 
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where V is a constant matrix of simple structure, e.g., diagonal, which is 
boundary compatible with K and L. Then for the boundary compatible set J = {V,K,L} 
containing the simple V matrix, it is often possible to analytically calculate 
the Green's matrices. We now have 


T J (y) = [I-U^rVoO-U^ y] 

where 

^ (y) = H J (t)C + yV(t,s){F(y(s),s) - Vy (s) }ds 
1 0 

U ^ L y = f G j (t,s){w(s) y (s) - Vy (s) }ds 
or 0 

1 

uj t y= J G J (t,s)6V(s)ds, 

0 

and finally we note 


3.5.29 


3.5.30 


3.5.31 


3.5.32 


T J (y) - y = H J (t)c + J G J (t,s){F(y(s),s) - WCs)y(s)}ds 3.5.33 

0 

so that j 

[ ( T y J ) ' - U^ L ] u = f G J (t,s){(3F/3y)(y(s),s) - W(s)}u(s)ds. 3.5.34 

0 

Hence we obtain the convergence benefits of choosing a general matrix W(t) while 
being able to calculate the Green’s matrices using the V matrix of simple 
structure . 

3.6. Applications of Contraction Mappings 

In this section we shall investigate the application of the contraction 
mappings method to the iterative solution of the TPBVP's arising from the regula- 
tion of nonlinear systems. In particular, using Theorem 3.4.14 we shall present 
the general form of the translated convergence theorems for the iterative 
solution of these TPBVP's. 
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Let us first consider the application to the system presented in Example 
2.3.1. Recall that this nonlinear system contained a nonlinear form containing 
only the state variable. We have for this case the following translated con- 
vergence theorem. 

Theorem 3.6.1. 

Let Yq(') be an element of ([0,1], R* 3 ) and let S = Sfy^.r). Suppose that 

(i) J = {V (t) ,M,N} is a boundary compatible set, and (iii) there are real 
numbers n and a with n > 0 and 0 < a < 1 such that 


1 


1) i|T J Cy 0 D - yfl = 


x o 

♦/ 

c J ct,si r 

'A(s) -B(s]R -1 (s)B' (s) 

' x Q (s) 



0 

0 

[ 

-Q(s) -A'(s) 

P 0 (s) . 


3.6.2 


V(s) 

x o (sr 

+ 

<?(x 0 (s)) 

Ids - 

f — \ 

4-» 

O 

X 




-P 0 (S) - 


3i|i 

3x (x 0 (s))p 0 (s) 

1 

-Po Ct) . 


» 


2) sup {Jam'll | 


sug 

sup 

1 

/ 

G J (t,sT 

'A(s) -B(s)R~ 1 (s)B(s) 

yes 

Bull <_ 1 

J 

0 

1 

-Q(s) -A'(s) 


[fii\ 

\3x/(x(s)) 


[D(x(s) ,p(s)) 




u(s)ds 


3.6; 3 


where D(s(s),p(s)) = [D k (x (s) ,p (s) ] = ^ ) (x(s)) p (s) , 

j.i k 1 

3) J— n < r 3.6.4 

1-a — 
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Then the CM sequence {y^ CO) for the TPBVP based on and J converges uniformly 
to the unique solution y* in S and the rate of convergence is given by 

n 


1 y* C - ) - y n (-)|li || yi (.) - y 0 (.)|| 


Proof: Apply Theorem 3.4.14 to the TPBVP of Example 2.3.1. 

From this general theorem statement, the performance of the numerical 
algorithm is difficult to predict. However, in the sequel, we shall develop 
coarse estimates for the convergence criteria contained in Theorem 3.6.1. 

We shall now apply the contraction mappings convergence theorem to the 
operator equation corresponding to the regulation of a system containing a 
general formulation for the nonlinearity, i.e., the TPBVP presented in Example 
2.3.45. The nonlinearity contained in that TPBVP is given as 


f(yCt)) = 


<Kx(t),p(t)) + B(t) ?(x(t),p(t)) 
-Z' (x(t) ,p(t))p(t) 


3.6.6 


where we defined 


u = C(x(t) ,p(t)) 

<t>(x(t) ,p(t)) = iHx(t) ,S(x(t),p(t)) 

and 

Z(x(t) >p(t)) = (3ip/dx) (x(t),£(x(t),p(t)) . 

Before applying the CM theorem, we shall first calculate an expression for 
(8f/3y) (y (t)) where y is the composite 2n vector [x,p]. We have 
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(3f/3y) (y) = 


(3<j>/3xXx,p) + B(t) (3£/3x)(x,p) 
- 3/3x [Z'(x,p)p] 


(3<(>/3p) (x,p) + B(t) (3C/3p) (x,p) 

3.6.7 

- 3/3p [Z'(x,p)p] 


Now defining the matrix functions D(x,p) and W(x,p) to be composed of the elements 

n 

D i;j (x,p) = ^ (32 ki /3x^) (x,p)p k 3.6.8 

k=l 

and 

WijCx.p) = ( az ki / 3 Pj) ( x .P)P k « 3.6.9 

k=l 

we have (3f/3y)(x,p) given as 


(3f/3y)(x,p) = 


Z(x,p) + B(t)(3C/3x)(x,p) 
-D(x,p) 


B(t) (35/3p) (x,p) + (3ij>/3p(x,p)| 


-W(x,p) - Z'(x,p) J 

3.6.10 


Using equation (3.6.10) we have the following theorem. 

Theorem 3.6.11; 

Let y Q be an element of ^?([0,1], R^) and let S = S(y Q ,r). Suppose that 
(i) J = (V(t),M,N) is a boundary compatible set for which (2.3.46) is differen- 
tiable, and (ii) there are real numbers n and a with n >_ 0 and 0 <_ ct < 1 
such that 
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1 


» |l J (y„)-y 0 | ■ 

At) 

x o 

* J G J (t,s) r 

' A(s) 

0 ■ 


V s > 



p .. 

° [ 

-Q(s) 

-As) 


p 0 (s) - 


-V (s) 


V s ?! 

lPo C s; > 


'<Hx 0 (s)',p 0 (s)) + B(t) e(x 0 (s),p 0 (s)J 
-Z'(x 0 (s) jPo Cs)) 


"j ds 


3.6.12 


r x 0 (tJ 

P 0 (t) - 


n , 


21 “I lF/>' 


sup sup 
yes |u|i <1 


X 

/ 


(t,s) [[ 


• A(s) 0 
Q(s) -A' (s)J 


-V(s) 


3.6.13 


2(x(s) ,p(s))+B(s) (3?/3x) (x(s) ,p(s)) 

-D(x(s),p(s)) 

B(s) (3?/3p) (x(s) ,p (s) ) + (3<(i/3p) (x(s),p(s)Jj J u(s)ds 
-W(x(s) ,p(s)) - Z*(x(s),p(s}) 




< a 


3) ^ r, < r 

Then the CM sequence {y n (-)} for the TPBVP based on y^ and J converges uniformly 
to the unique solution y* in S and the rate of convergence is given by . 

y*0) - y n (-) I yjO) - y 0 (O • 3.6.14 

As we have indicated, cursory examination of Theorems 3.6.1, 3.6.10, and 
3.6.21 yields limited information conveming the convergence of the CM sequence. 
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The difficulty to a great extent lies in the intricacy of evaluating the integral 
containing the Green's function, G^t.s), and the derivative term of the form 
(3F/3y) (y (s) ) - V(s). In the next chapter, we shall consider techniques for 
alleviating these difficulties so that meaningful convergence analysis can be 
made without extensive computation. 


, ~ s 
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CHAPTER 4 


CALCULATION OF CONVERGENCE CRITERIA 

4.1. Introduction 

For the boundary compatible set J = (V(t) ,M,N} , we consider the iterative 
solution of the operator equation 

X = T J (y) 4.1.1 

where T^(y) is given by 

y(t) = T J (y)(t) = H J (t){c-g(yCO) - h(y(l)) + My(0) + Ny(l)} 

! 4 - 1 ' 2 
+ J G J (t,s){F(y(s),s) - V (s)y (s) }ds 
0 

and the Green's functions H* 1 , G J are given by 
H J (t) = 4> V (t,0)[M + N4 V (1,0)] -1 

and 

f G I J (t,s) = 4 V (t,0)[M + OT> V (l,0)]~ 1 M«> V C0,s) , 0 <^s < t 

4.1.3 

G^(t,s) = -<J V (t,0)[M N4> V (l,0)]” 1 N4> V (l,s) t < s <_ 1 . 

Theorem 3.4.14 specified conditions necessary for convergence of the CM sequence 
y n+ j = T^(y n ). In this chapter, we discuss in detail the evaluation of the 
convergence criteria. In particular, we discuss two general schemes that may be 
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used to lessen the analytical difficulties involved in calculating the convergence 
parameters n and a. 

The first scheme is simply that of selecting- very simple V matrices for use 
in the representation. For example, one might select V as the zero matrix or a 
constant diagonal matrix. For these matrices the fundamental matrix is readily 
obtained and the Green's function matrices are often easily calculated. 

The second scheme involves the use of a similarity transformation. In this 
approach, a more general constant V matrix is selected and transformed into a 
canonical fonn. Then using the canonical form, the fundamental matrix is obtained. 
However, for this approach, the calculation of the Green's function matrices is 
somewhat complicated by the transformation matrices. In conclusion, an approximate 
technique is developed which often yields accurate estimates. 

4.2. Estimates of Convergence Criteria 

Before considering specific boundary compatible sets, we first specify those 
estimates of the convergence parameters which are desired. As indicated in 
Theorem 3.4.14, the numbers to be calculated are estimates for ||t J Cy Q D ~yQ || 
and |(T y J r||. 

First consider the estimation of |t J C y Q D _y q II ' At t * lis P°i nt > it will be 
useful to discuss an effective technique for obtaining the initial estimate of 
the solution. Consider the iterative solution of the nonlinear TPBVP 

y = F(y,t) 4.2.1 

Ky(0) + Ly(l) = c, 

and the choice of the boundary compatible set J = {W(t),M,N} to be used in the 
integral representation. Since the boundary conditions of (4.2.1) are linear. 
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we choose M=K, N=L in the representation. If we now choose W(t) based upon 
(3F/3y) (y,t) , i.e., a linearization of the system, then the solution to the 
linear TPBVP 


f = W(t)y 
Ky(0) + LyCl) - c 


4.2.2 


is often a good initial estimate for the solution of (4.2.1). Moreover, this 
choice considerably simplifies the calculation of T J (y 0 )-y Q since y 0 (t)=H J (t)c 
and 1 


r Vo>-^o ■ / G J (t,s){F(y 0 (s) ,s) - W(s)y Q (s) }ds 


4.2.3 


0 

for the boundary compatible set J = {W(t),M,N}. 

The other norm which must be calculated is the derivative norm ||(T ^)'| • 
As presented previously in (3.3.7), a coarse estimate for ||(T^ J ) ' || is given as 




sup sup | 
i£P t ( 

j=l 0 

M)\ 



, P 

E 

l k=l 


4.2.4 


Let us make the following definitions. 

Definition 4.2.5. 

Let P(t) = [p^j(t)] be a matrix with entries 


■ x 

p i j C t ) = f \g\j (t,s) |ds 
0 


4.2.6 


P.-Ct) = f |g J T (t ,s) | ds + f \gt (t,s) | ds 

J J 0 % { n ij 


4.2.7 
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where g^ and are elements of G^ft.s) and Gjj J (t,s) as given in 

ij ij 

Definition 4.2.8 

Let z = [z ] be a vector with elements 
i P 


z 0 = SUp 
U i t£ [0,1] < 1 u 


£ v ij (t) v (t) i| 

j=l 3 


(4.1.3). 


4.2.9 


Definition 4.2.10 

Let z = [ z . ] be a vector with elements 
L 

/ P ) 

z = sup sup £| (3F /3y ) (y(t),t) - v. . (t) | . 4.2.11 

1 t£[0,l]y£S{^ 1 3 13 ) 

From (4.2.3) and (4.2.4) it follows that conservative values for the convergence 

parameters n and a are given by 


and 




4.2.12 


4.2.13 


In the remainder of this chapter we shall be primarily concerned with techniques 
for determining the matrix P(t) for boundary compatible sets containing simple 
V matrices. 


4.3. Boundary Value Sets of Interest 

In this section we shall briefly specify the form of those pairs of boundary 
condition matrices which are of interest. The necessary conditions for regulation 
of nonlinear systems reduced to TPBVP's of the form 


y = Sy + f (y) 4.3.1 

My (0) + Ny (1) = c 4.3.2 
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where the matrices M and N depended on the quadratic cost functional being used. 
Specifically we had the following cases. 

Definition 4.3.3. 

For quadratic cost functionals including a terminal state penalty of the 
form (x(T),Kx(T)) , the boundary condition matrices were 


'I 

0 

N = 

‘ 0 

O' 

p 

0 


-K 

I. 


Since we have rank [M N] = 2n, Lemma 3.2.8 assures a matrix V exists so the set 
J = {V,M,N} is boundary compatible. We shall henceforth refer to set (4.3.4) as 
boundary value set {1}. 

Example 4.3.5. 

For quadratic cost functionals which do not include a terminal penalty, 
the boundary condition matrices are given as 


I 

o' 

[0 

o' 


M = 


N = 



0 

0 

Lo 

i. 

4.3.6 


Again rank [M N] = 2n, so a V matrix exists such that J = {V,M,N} is boundary 
compatible. The set (4.3.6) shall be referred to as boundary value set (2). 

4.4. Boundary Set for Regulation with Terminal Cost 

In this section the use of simple V matrices with boundary value set {1} 
will be considered. The requirements for boundary compatibility of the various 
sets J = (V,M,N} will be noted in particular. 

Boundary value set {1} is given specifically as 
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4.4.1 



'I O' 


o 

o, 

M = 


N = 



o' 

o 


--K I. 


and a general 2n x2n V matrix is represented as 


V = 


V 

V 


11 

V T2 

21 

V 22 


4.4.2 


The fundamental matrix for V is represented as 


V 

* (t,s) 


fi 11 (t,s) fl 12 (t,s)' 

L« 21 (t,s) fl ? ,(t,s) 


22 ' 


4.4.3 


The matrix [M + N*^(1,0)] is now formed explicitly as 

I 0 


M+N<t> (1,0) = 


l-Kn n (l,0) + J2 21 Cl,0) -Kfl 12 (l,0)+fl 22 Cl,0) 


4.4.4 


and the inverse, if it exists, may be written as 


[M+N$ V (1 ,0) ] ~ 1 = 


I 

-[-Kft 12 (i,o) + n 22 (i,o)]" 1 [-Kfi 11 (i,o)+n 21 (i)]' 

0 

[-K£i 12 (i,o)+a 22 ci,0)] _1 _ 

4.4.5 


For this inverse to exist, the matrix [-Kfi 12 (l,0)+& 22 (1,0)] must be nonsingular. 

It is noted that for V equal to the zero matrix or a diagonal matrix, the set 
J = {V,M,N} is boundary compatible. The core of the Green's function is given 
by the matrices [M+N4 V C1,0)] and [M+N$ V (1,0)] which are explicitly given as 
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[M+N$ V (1,0)]" 1 M = 


I 0 

,-l 


L-[-Kfl 12 U ,0)+a 22 ci,0)] [-Ko n (i,0)*tt 21 )i,0)j OJ 

4.4.6 


and 


[M+N4 (1,0)] 




.-1 


Kfl 12 (i,0)*n 22 (i,0)] k [-icn 12 )i,0)+n 22 (i 


.o)]' 1 ] 


4.4.7 


We shall now consider specific choices for the V matrix. 

Example 4.4.8. 

Consider the choice of the simplest V matrix, i.e., assume V = 0. The 
fundamental matrix is then given as 


V 

* (t,s) = 


I 0 
0 I 


4.4.9 


Now using (4.4.6) and (4.4.7), 


and 


[M+N$ V (1,0)] _1 M 


[M+N4> V (1,0] _1 N = 


I 0 

K 0 . 

0 O' 
-K I 


4.4.10 


4.4.11 


The Green's function matrices are calculated as 


Gj(t,s) = $ V (t,0)[M+N« V (l,0)] -1 M$ V (0,s) = 


I 0 
K 0 


and 


(t,s) = -« V (t,0)[M+N* V (l,0)] _I N* V (l,s) = 


0 0 
K -I 


4.4.12 


4.4.13 
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and the 2n x 2n P(t) matrix defined as 


L A. 

P(t) = / | G j (t , s ) | ds + /|G^(t,s)|ds 


4.4.14 


is calculated to be 


P(t) = 


tl 0 

I K I Cl-t)I. 


4.4.15 


where elements of | K J are given as |k_|. 

Example 4.4.16 

The use of a p x p (2n x 2n) diagonal V matrix is now considered. Let V 
be represented as 


V = 


n J 


4.4.16 


so the fundamental matrix is then simply 


* (t,s) = 


X 1 (t-s) 

e 

x n (t ’ s) 

e 

— 

, 

0 


P, (t-s) 

e 1 

• 

• 

0 

e 


4.4.17 
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which shall be denoted as 


$ (t,s) = 


fi n (t,s) 0 


0 S 22 (t,s) 


This yields using (4.4.6) and (4.4.7), 


[M+N4> V (1 ,0) ] -1 M = 


fi 22 (l,0)K Q n (l,0) 0. 


and 


[M+N4 V (1 ,0) ] ” 1 N = 


-£2^(1,0)K «'*(!, 0) 


The Green's function matrices are determined to be 


G j (t , s ) = 


Q n (t,s) 




Lfl 22 (t,0)n 22 (i,0)KR n (i,s) 0 


and 


Gj j (t , s) = 


0 


,-l 


L o 22 (t,o)n 2 *(i,p)Ko n (i,s) -n 22 (t,s)J 


In many instances the K matrix associated with the terminal cost is a 
diagonal matrix. Let us now assume K diagonal with elements k^. Then using 
(4.4.21) and (4.4.22), the P(t) matrix is found to be 


4.18 


4.19 


4.20 


4.21 


4.22 
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P(t) = 


1 V 
T O -W 

> N 


\ 


■' 1 f V n 
- (e - 1} 
n 


A i 

- 7 — e ( 1 -e ) 


\ 


k -y (1-t) X 
1 n 1 H n ,, n N 
1 — e (1-e ) 


n 


i n -"it 1 -*) 

— (1-e 
y l 

\ 


) 


\ 


\ 


1 n 
— (1-e 

% 


4.4.23 


Example 4.4.24. 

Many nonlinear systems of interest have an underlying oscillator structure. 

For this reason we shall consider a choice of V matrix containing linear oscillator 
elements. This choice is represented as 

S, 

• r 

0 


V = 


S. 

J 


4.4.24 


where the are 2x2 matrices of the form 


S. 

l 


a. to. 

l l 


-o). a. . 

i i J 


4.4.25 


The fundamental matrix for this choice of V is given as 
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or 


(t,s) = 


* (t.s) = 


♦j(t.s) 

• 


• 

* 4> j (t,s) 

0 

-- 

<f> k (t,s) 

• 

1 

O 

• 

‘ ♦„<*-*> . 

ra n (t,s) o 



fi 22 Ct»s). 


4.4.26 


4.4.27 


where the 4 >^(t,s) are 2x2 matrices of the form 


^(t.s) = 


a i (t-s) 

e cos uK(t-s) 

a i (t-s) 

-e sin (ik (t-s) 


a i (t-s) 

e sin (t-s) 


a . 


* *- t_s ^cos u. (t-s) 


4.4.28 


From (4.4.5), the matrix [M+N<t V (1 ,0) ] is nonsingular if ft 22 (l,0) is nonsingular. 
We now have 


where 


fi 22 (l,0)= 


.-1 


( 1 , 0 ) 


Vd’°). 


^ 1 C1,0)'- 


-a . 

e 1 cos a). 

i 

-a. 

e 1 sin to. 

l 



sin to. 
i 

cos to. 
1 


so this choice leads to a boundary compatible set. 
The Green* s functions are found to be given as 


4.4.29 


4.4.30 
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Gj (t,s) = 


D n (t,s) O' 

Lfi 22 (t,0) A Q n (0|S) 0 


4.4.31 


and 


Gj T Ct,s) = 


4.4.32 


1 1 

l 

-n 22 (t,o) a fi n (o,s) -n 22 (t,s)J 

> 


where the matrix A 

is given as 



a = -« 22 (i,o; 

The matrix P(t) is thei 

) K 0 U (1,0). 

i given as 
t 


4.4.33 

P(t) = 

J 

f|n i]L (t,s) |ds 
' 1 

fr r 

0 


c 

1 

D 22 (t,s) | ds 

4.4.34 


■ J 

c 

j |n 22 (t,o) A D n (0,s)|ds J | 

1 0 



Due to the oscillatory nature of the elements of G (t,s), the integration of 
the absolute values somewhat complicates an analytic solution for P(t). However, 
in a future section we shall consider approximate techniques for obtaining this 
P(t) matrix. 


4.5. Boundary Set for Regulation with No Terminal Cost 

In this section the use of simple V matrices with boundary value set {2} 
is considered. The requirements for boundary compatibility of the various sets 
J = (V,M,N,} shall be noted in particular. Boundary value set {2} is given 
specifically as 


'i 

O' 


'0 

o' 



N = 



.0 

0 , 


0 

I 
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A general 2n x 2n matrix is represented as 


V = 


V V 

11 12 

V V 

21 22 


4.5.2 


and the corresponding fundamental matrix is given as 


V 

* (t,s) = 


« 11 (t,s) n 12 (t,s) 

Ln„(t,s) fi„(t,s)J 


4.5.3 


The matrix [M+N$ V (1,0)] is formed as 

I 


[M+N4 (1,0)] 


Ln 21 Ci,o)- fi 22 (i,o) 


4.5.4 


and the inverse, if it exists, is given by 


[M+N$ V (1 ,0) ] _1 


L -«‘^(l,0)ft 21 (l,0) 


»-]u) 


4.5.5 


For this inverse to exist, il 22 (l,0) must be nonsingular. It is noted that for 
V equal to the zero or diagonal matrix, the set J = (V,M,N) is boundary compatible. 
At this point we shall begin to take advantage of the fact that the remaining 
results desired in this section may be obtained from the results of the previous 
section with K equal to zero. These results are now presented for V matrices 
of simple structure. 

Example 4.5.6. 

The first selection for the V matrix is the zero matrix, i.e., V = 0. 

Using the results of Example 4.4.8 with K = 0, the P(t) matrix is given as 
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P(t) 


tl 0 

0 (l-t)I 


4.5.6 


Example 4.5.7 

For the second choice of the V matrix, a p * p (2n x2n) diagonal matrix is 
selected, i.e., 

N , 

0 


V = 


4.5.8 


Now specializing the results of Example 4.4.16 with K =0, the P(t) matrix is 
obtained as 

" i V 
~ (e -1) 

1 


P(t) = 


1 r * nt l 

r- (e -l) 
n 


1 n 

(1-e ) 


y l 


1 n 
— (1-e 

y n 


4.5.9 


4.6. Application of Similarity Transformations 

As an introduction to the use of similarity transformations, consider the 
linear TPBVP 


y = Vy My (0) + Ny(l) = c . 


4.6.1 
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If the set J = (V,M,N) is boundary compatible, the solution to (4.6.1) is given 
by Theorem 3.23 as 

y (t) = 4> V (t,0)[M+N4> V (l,0)]' 1 c. 4.6.2 

V 

In an attempt to ease the calculation of the fundamental matrix 4> (t,0), consider 
the use of the nonsingular linear transformation 

A.z = y . 4.6.3 

From (4.6.1) , the transformed TPBVP is given as 

i = A _1 VA z, MAz(O) + NAz(I) = c . _ 4.6.4 

If the set J = {A ^VA > MA, NA) is boundary compatible, the solution for (4.6.4) 
may be written as 

z(t) = $ A VA (t) [MA + NA$ A VA (1,0) ] -1 c. 4.6.5 

In passing, it may be quickly shown that if the set J = (V,M,N) is boundary 

compatible, the transformed set J = {A* VA, MA,NA) is also boundary compatible. 

V 

With the matrix [M+N4> (1,0)] nonsingular, post multiplication by A yields the 

nonsingular matrix [MA+N$ (1,0)A]. Th,e fundamental matrices are related by 

V A'^VA, -1 V 

$ (1,0) = A4> (1,0)A so the nonsingular matrix [MA+N$ (1,C)A] may be written 

as [MA+Nt |A ^ A (l,0)] indicating the transformed set J = (A *VA,MA,NA) is boundary 

compatible. If the transformation A *VA reduces V to a canonical form, the 

fundamental matrix $ A VA (t,s) is of simple structure. 

Now consider the nonlinear TPBVP 

y = Sy + f(y) My(0) + Ny(l) = c . 4.6.6 
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Again consider the nonsingular linear transformation 


A z = y 


4.6.7 


and let 


D = A" L VA . 


4.6.8 


Then (4.6.6) becomes the transformed TPBVP 


.-1 -1 
z = A SAz + A f(Az) 


4.6.9 


MAz(O) + NAz(l) = c . 


•v r *“ 1 

If the set J = {A VA, MA, NA} is boundary compatible, the integral representation 
for (4.6.9) is 

1 

•v <v _ V 

T J (y) = tf^^c + J G J (t,s){A' 1 SAz + A' 1 f(Az(s)) - Dz)ds 4.6.10 

0 

where the Green's functions are given as 

H J (t) = $ D (t,0)[MA + NA$ D (l,0)}' 1 4.6.H 


and 


T i 4 D (t,0)[MA + NA3> D (l,0)] _1 M4> D (0,s) , 0 < s < t 

G Ct ’ s3 = ' D DID " " 

(t,0) [MA + NA$ U (1 ,0) ] W'Q.s) , t < s < 1 . 


4.6.12 


If it is desired to investigate the iterative solution of the operator equation 


z = i (z) , 


4.6.13 


the operator derivative (T ) 1 is given as 

1 

(T z J )'u = / G J (t,s){A _1 SAz (s) + A _1 (3f/3y)(Az(s))A - D}u(s)ds 4.6.14 

0 
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if the TPBVP is differentiable. However, rather than using (4.6.14), another 
approach may be taken. It may easily be shown that a direct relationship exists 
between the Green's functions for the boundary compatible set J = (V,M,N) and 
the transformed boundary compatible set J = (D,MA,NA) . In particular, 

H J (t,s) = AH J (t,s) 4.6.15 

and 

G ^ (t , s ) = AGj(t,s)A‘ 1 4.6.16 

G^(t,s)= AG J nCt,s)A^ 4.6.17 

Hence the integral representation for (4.6.6) may be written as 

1 

y(t) = T* 1 (y) (t) = AH J (t,s)c + y AG J (t ,s) A _1 {Sy (s)+f (y (s) ,s) -Vy (s) }ds 

0 4.6.18 

Then if the matrix A _1 VA is a canonical form, $ D (t,s) and G J (t,s) are often much 
easier to calculate, and it may very well be easier to calculate estimates for 
(T J ) ' . 

v y 

The theory of canonical forms has received great attention in the past years. 
General books of interest include Gantmacher [Gl], Bodweig [B2] .Turnbull [Tl], 
and Ferrar [F2], Of interest to control analysts are the books of Bellman [Bl] 
and Ogata [01]. In particular, we now present a well known theorem concerning 
the diagonalization of matrices. 

Theorem 4.6.19 

If the characteristic roots of the matrix V are distinct, there exists 
a matrix A such that 
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4.6.20 


Proof. (See Bellman or Ogata) . 

However, if a p * p matrix V does not possess p linearly independent eigenvectors, 
then V is not similar to a diagonal matrix. In this case, it can be proved 
rigorously that a p x p matrix, V, possessing less than p linearly independent 
characteristic vectors is similar to the Jordan canonical form, where the 
elements in the main diagonal are the characteristic roots and the elements 
immediately above the main diagonal are either one or zero and all other elements 
are zero. (The proof of this statement may be found in Turnbull.) However, 
rather than using the more involved Jordan canonical form, we shall make use of 
the following result from Bellman. 


Theorem 4.6.21 

Given any matrix W, we can find a matrix V with distinct characteristic 
roots such that ||W-Vll C , where C is any preassigned quantity. 

Proof. (See Bellman.) 

The importance of Theorem 4.6.21 is as follows. Assume analysis of the 
convergence conditions indicates the matrix W is a good choice for use in the 
integral representation. If W contains multiple characteristic roots, it is 
not similar to a diagonal form and the advantages of this simple form are not 
available. However, since we are free to choose the matrix, we may use Theorem 
4.6.21 and "perturb" the W matrix to a V matrix "close to W" (i.e., ||W-V|| <_ C ) 
which does have distinct characteristic roots. We may then determine a matrix 
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A such that A ^VA is a diagonal form. We shall consider one further special 
case and that is the system whose distinct characteristic roots involve complex 
conjugate paris. We shall follow Ogata. 

Assume for convenience that the system involves only one pair of complex 
conjugate characteristic roots. Extension to the case where there are more 
than one pair of complex conjugate characteristic roots is obvious. Assume 
that the eigenvalues X^ and X^ are complex conjugates and sire given by 

Xj = ct + j o> Xj = a - ju . 4.6.22 

Assume also that the eigenvalues X_,...,X are real and distinct. The diagonal 


4.6.23 


4.6.24 


4.6.25 


4.6.26 


matrix is then of the form 

a+jw 


D = A' 1 VA 


CJ- JO) 


o y 

which may be transformed into the modified diagonal form 


D = 


a w 
-u a 


3 ~ - 


0 


V 


by means of the transformation matrix defined by 

ol 


K = 


1/2 - j/2 

1/2 j/2 


.1 


Namely the modified diagonal form D is given as 


D = K" 1 DK = K -1 A' 1 VAK . 
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Not only does D have only real elements but, more significantly, and 

AK have only real elements. Now setting 

A = AK 4.6.27 

?-l „-l.-l 

A = K A , 

the transformation is given the standard form of 

D = A' 1 VA . ' 4.6.28 


As a result of this discussion, in the sequel the term canonical form shall 
specifically refer either to diagonal or to the modified diagonal form as in 
(4.6.24) . 

We shall now choose several forms for the V matrix to illustrate the use 
of the similarity transformation with the boundary value sets of interest. 
Example 4.6.29. 

Let us consider the boundary value set 


I 

0 


■ 0 

O' 



N = 



.0 

0 


-K 

I 


and the 2n x 2n V matrix with distinct characteristic roots 


V = 


V 11 ° 


22 


The similarity transformation has the form 


A,, 

0 


r-H 

1 

< 

0 

11 

.-1 

11 



A 


-1 

.0 

a 22 - 


.0 

A 22 


4.6.31 


4.6.32 
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and the canonical matrix D is given as 
D = A _1 VA 


The fundamental matrix of the canonical matrix has the general form 


* D (t,s) 


o ■ 

0 <)> 22 Ct,s) 


The matrix [MA + NA4>^(1,0)] is obtained as 

A, 


[MA + NA$°C1,0)J = 


11 


-KA^^M) A 22 * 22 (1,0). 


and if the inverse exists. 


[MA + NA* D (1,0)] _1 


l ll 


♦ 2 J(1,0)A 22 KA n ^ n (1,0) >"Jci,0)A 2 J J 


The Green's functions are then found as 


Gj(t,s) 


4' 11 (t,s) O' 

< t’ 22 (t ,0) A<}> 11 (0,s) 0- 


and 



0 

. l i> 22 (t,0)A(fi 11 (0,s) 


0 

- < (' 22 (t,s). 


where the n x n matrix A is given by 


4.6.33 


4.6.34 


4.6.35 


4.6.36 


4.6.37 


4.6.38 
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4.6.39 



4.6.40 


4.6.41 


Example 4 S 6.42. 

Let iis now consider the use of a 2n * 2n V matrix of the form 


4.6.42 


with the boundary value set corresponding to no terminal penalty, i.e.. 




4.6.43 


Defining the canonical form D = A~*VA, we can calculate P(t) by specializing 
the results of Example 4.6.29 with K = 0. P(t) is obtained as 
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In this, section, the use of similarity transformations was introduced in an 
attempt to simplify the calculation of P(t). For some cases the technique worked 
very well yielding simple expressions for P(t). However, in some instances the 
matrix P(t) is very awkward to calculate. Consideration of the similarity trans- 
formation led to the development of an approximate technique which is presented 
in the next section. 
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4.7. Approximate Technique 

We shall now introduce a technique which, though not mathematically rigorous, 
allows one to obtain estimates for P(t) in a much simpler fashion. Using the 
canonical transformation A, define the matrix D as D = A _1 VA and write G^(t,s) 
and Gjj(t,s) as 

Gj(t,s) = {A4 D (t,0)A" 1 H[M+N$ V (l,0)]' 1 M}{A4 D (0,s)A' 1 } 4.7.1 

Gjj (t,s) = -{A4 D (t,0)A' 1 }{ [M+N<J V (l,0)]‘W / (l,0)}{A'f D (0,s)A' 1 }. 4.7.2 


The terms have been separated by brackets to indicate the factors contributing to 
the magnitude of P(t), namely the inversion and the integration of the fundamental 
matrices. Consider the general 2n x 2n V matrix and the resultant fundamental 
matrix to be given as 



V 11 

V 12 _ 

v 

^ n (t,s) 

fl 12 (t,s) 

V = 

V 21 

CN 

CM 

> 

$ (t,s) = 

_n 21 (t,s) 

n 22 (t,s) 


For the boundary value matrices 


i 

o' 


" 0 

O' 



N = 



.0 

0. 


,-K 

I. 


we have 


[M+N4 V (1,0)] 


I 0 

_-Kfi n (i,0) + n 21 (i,0) '. K n 12 (i,0)+n 22 (i,0) 


4.7.4 


4.7.5 


and if the inverse exists. 
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[M+N^Cl.O )]" 1 


I 

. -[-Kfi 12 (i ) 0)+n 22 (i,0)]" 1 [-Kn n (i. 0 )+n 21 (i, 0 )j 


o 

[-icn 12 (i,o)+n 22 (i,o)] 


4.7.6 


The center bracketed terms are then found as 


[M+N« V (1,0)] _1 M = 


I 

A 


0 

0 


and 


[M+N$ V (1 ,0) ] 1 N$ V (1,0) 


0 0 
-A I_ 


where the n * n matrix A is given as 


4.7.7 


4.7.8 


A = -[Kn 12 (i,0) + n 22 (i,0)]" 1 [-Kn 11 (i,0)+n 21 (i,0)]. 4.7.9 

Now assuming that $ D (t,s) represents the primary magnitude characteristics of 
A4> D (t,s)A 1 , we shall form 


{4> D (t,0) }{ [M+N$ V (1,0)] _1 MH$ D (0,s)} 4.7.10 

and 

-{<f D (t,0)}{ [M+N4> V (1 ,0) ] _1 N$ V (l,0)}{4 D (0,s)} 4.7.11 

as approximations to Gj(t,s) and G^(t,s). Following the discussion in Section 4.6 
concerning canonical forms, the fundamental matrix <J^(t,s) may be represented in 
the form 
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$ D (t,s) = 


K 11 


(t,s) 


0 


22 


(t,s) 


4.7.12 


Using (4.7.10) and (4.7.11), we obtain the approximations 

ol 


Gj(t,s) 


^jCt.s) 


^ 2 Ct,0)&4,J 1 (0,s) 0 


4.7.13 


and 


l! (t ,s ) * 


^22 ^t»0)A<ji^ 1 (0,s) 


Finally, this yields 




P(t) 


U?, (t,s)|ds 


/' 

J'l I ^22 C 1 >°) A *11 c°> S D I ds M 2 «,sU*s 


4.7.14 


4.7.15 


For the boundary value set 



I 

0 


0 

O' 

M = 



. N = 




0 

0 


.0 

I 


4.7.16 


we may specialize the results of (4.7.15) with K = O' to obtain 
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Q 


P(t) 



1 

J l dS 

*t 


4 . 7.17 


where 

1 = • 4 . 17.18 

These approximations greatly simplify the calculation of P(t), and moreover, they 
capture the primary quantitative behavior of the P(t) matrix. The concepts and 
techniques introduced in this chapter will be illustrated in several numerical 
examples in Chapter 6. 
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CHAPTER 5 


CONTROLLABILITY FOR NONLINEAR SYSTEMS 


5.1. Introduction 

The concept of null controllability is a natural aspect of the study of 
optimal regulation for nonlinear systems. Whereas the optimal regulator attempts 
to drive the system from its initial state into a region near the origin, null 
controllability is concerned with driving the system precisely to the origin. 
Historically, the issues of regulation and controllability are closely intertwined. 
The study of linear regulator problems in a general framework served to uncover 
some of the underlying relationships that exist between the structure of the 
optimal system and the fundamental concept of controllability [Kl] , [K3] . Much 
of the effort to date concerning null controllability of nonlinear systems has 
involved determination of feedback controllers such that the driven systems 
satisfy certain Lyapunov-type stability arguments [B4], [G3] . In this chapter 
the integral representation of TPBVP's and the contraction mapping theorem will 
be used to investigate the controllability of nonlinear systems via existence of 
solutions arguments . 

5.2. Controllability for Linear Systems 

As an introduction to the controllability issue and the approach to be 
taken in the study, we shall first consider the controllability of linear systems. 
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Definition 5.2.1. 


The autonomous linear control process 


x(t) = Ax(t) + Bu(t) . 5.2.2 

with u£fi = R m , is (completely) controllable in case: for each pair of points 

Xg and Xj in R n , there exists a bounded measurable controller u(t) on some finite 
interval 0 <_ t <_ T , which steers Xq to x^ . 

Theorem 5.2.3. 

The autonomous linear system 

x(t) = Ax(t) + Bu(t) 5.2.4 


with u£i! = R m , is (completely) controllable if and only if a solution exists to 
the linear TPBVP 



Proof: Assume a solution x*(t), p*(t) exists to the TPBVP (5.2.5). Now consider 
the optimization problem of determining a control u(t) to drive the system (5.2.4) 
from the initial state x^ to the terminal state x^ such that the cost functional 



5.2.6 
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is minimized. Application of the Pontryagin principle yields precisely the 
TPBVP (5.2.5). Then the control 


u(t) =-B'p*(t) 5.2.7 

drives the system from x Q to x^ . 

Conversely, assume that the system (5.2.4) is completely controllable. We 
shall show that a solution exists to the TPBVP. The linear TPBVP 


y(t) = Vy (t) + f(t) 5.2.8 

My (0) + Ny (1) = c 

with 



"a -bb 1 ' 


I O' 


0 O' 

V = 


M = 


N = 



.0 -A' 


o' 

o 


.1 0. 


V 

has a solution for every f(t) and c if det[M+N4> (1,0)] f 0 . The fundamental 
matrix for V is given in the form 


;n (t,s) 0 (t,s)-, 

4> (t,s) = 1 11 1 

n 22 ( t » s ) 


n ll (t 

0 


5.2.10 


and the matrix [M+N$^(l ,0) ] is obtained in the form 


[M+N4 (1,0)] = 


n n U,o) n 12 (i 


,o)] . 


5.2.11 


The inverse, if it exists, is given as 

I 


[M+N<J> V (1 ,0) ] _1 


-n' 2 (i,o)fi n (i,o) 


5.2.12 
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and det[M+N4^ (1 ,0) ] = det [fl 2 (1,0)] . Now investigating the differential equation 

describing 4^(t,0), we have 


si Ct,0) = An ,(t,0) - BB'n (t,o) , n..(o,o) = o 

1 1 w 5.2.13 

o 22 (t,0) = -A*n 22 (t ,0) , fi 22 (0,0) = I 

These equations yield 

n 22 (t,0) = 4>" A ' (t,0) = $ A (0,t)' 5.2.14 

and 

n i2 (t,0) = <t>' A (t,0) //(O^oiBB./fO^'do . 5.2.15 

0 


Hence for the existence of a solution to the TPBVP (5.2.5), we must have 

1 

det[4> A (l,0) f $ A (0,a)BB , 'I> A (0 > o) 'da] /< 0 . 


5.2.16 


However, the assumption of complete controllability specifies 
1 

det [ y $ A (O,0)BB'4> A (O,o) ’do] { 0 , 

0 


5.2.17 


therefore a solution to the TPBVP must exist. 

Hence the issue of invertibility of n^ 2 (l,0) leads to the well known 
controllability Grammian and the approach is seen to yield conditions compatible 
with previously derived results. The following corollary will often prove useful 
when selecting a V matrix for a controllability investigation. 


74 



Corollary S.2.18 


Let the constant 2n x 2n matrices V,M, and N be of the form 


V 11 

V 12 

M = 

I 

O' 

N = 

0 

O' 

0 

- V lJ 


0 

0 


I 

0. 


with = -BB' where B is an n x r matrix. Then the set J = (V,M,N) is 
boundary compatible if and only if 

rank [B,V n B, . . . .v"" 1 B] = n 5.2.20 

Proof. See proof of theorem 5.2.3 and Lee and Markus [LI] for the relationship 
between (5.2.20) and the controllability Grammian. 

The obvious advantage provided by Corollary 5.2.18 is that it removes the 
calculation of $ V (1,0) when determining the boundary compatibility of a set J 
in the form of (5.2.19). With this background, we shall now consider nonlinear 
controllability. 

5.3. Nonlinear Controllability 

In this section we shall extend the approach of Section 5.2 to include 
nonlinear systems. However, rather than considering global controllability 
as for linear systems, we shall consider local null controllability, i.e., 
the problem of regulating an initial state, near the origin, to the origin. 
Definition 5.3.1. ['Ll] 

Consider the control process in R n 

x = f(x,u) in in R n x ft 5.3.2 

where ft is a restraint set in R m . The domain 3f null controllability is 


75 



0 


defined as the set of all points x^ £ R n , each of which can be steered to = 
by some bounded measurable controller u(t) CO in finite time. If e/f^contains 
an open neighborhood of x^ = 0, then (5.3.2) is said to be locally controllable 
(near the origin). We shall now consider the null controllability of nonlinear 
systems by means of integral representations . 


Theorem 5.3.3. 

Consider the control process in R n 

x = f(x,u) in C 2 in R n * 0 

with u = 0 interior to the restraint set 0 C R m . 

Assume 

(a) f(0,0) = 0 

(b) rank [B,AB,...,A n 1 B] = n 

where A = (3f/3x)(0,0) and B = (3f/3u)(0,0) 

Then the domain*^ 3 af null controllability is open in R n . 
Proof. Let us define the function \Ji(x,u) as 

iji(x,u) = f(x,u) - Ax - Bu 

so that 

¥( 0 , 0 ) = 0 

(J|) (0,0) = 0 

and 

(|i)(0.0) - o . 

Now consider the optimization problem composed of the system 


5.3.4 


5.3.5 


5.3.6 


5.3.7 


5.3.8 


5.3.9 


5.3.10 


5.3.11 
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x = Ax + Bu + i|/(x,u) , 


5.3.12 


the boundary conditions , 

x (0) = x 0 , x(lj = 0 , 

and the cost functional 
1 

J = j f <u(t),u(t)) dt. 
0 


5.3.13 


5.3.14 


The Hamiltonian for this problem is given by 

H = j <u(t),u(t)> + (Ax(t),p(t)) + (Bu(t).pCt)^ +' (<MX,u),p(t)) 

5.3.15 

and the costate variable is described by the differential equation 

P = A ’P - (J|)’ Cx.ujp . 5.3.16 

Now assume a control of the form 

u = -B'p 5.3.17 


accomplishes the desired transfer. The canonical system of equations is now 
given as 

x = Ax - BB'p + 1 J 1 (x ,u (p) ) 5.3.18 

£ = -A'p - (||)' (x,u(p))p 5.3.19 

subject to the boundary conditions 
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5.3.20 


I o' 


x(oy 

+ 

o 

o 

xur 


o 

X 

O 

O 


P(0) 


Li o J 

p(i). 


. o j 


This may be expressed as 


where 


and 


y = Sy + F (y ) 
My(0) + Ny(l) = c 



A 

-BB' 


I 

0 ' 


0 

0 


x o‘ 

s = 



M = 



N = 



C = 


0 

-A' . 


. 0 

0 


I 

0 . 


.0 . 


F(y) 


>Kx,u(p)) 
(H-) 1 (x,u(p))p 


5.3.21 


5.3.22 


5.3.23 


For the boundary compatible set J = (V,M,N) where M and N aTe given in (5.3.22), 
the solution to (5.3.20) may be written as 


1 

y(t) = tf^tjc + J G J (t,s)(Sy(s) + F (y (s ) ) - V(s)y(s)}ds 5.3.24 

0 

or as an operator equation 

y = T J (y). 5.3.25 

Clearly with c = 0, y(-) = 0 is a fixed point of T . We are now interested in 
the existence of a solution, y(-), if c is varied in a neighborhood of the origin. 
Define the variable z as 

z = P J (c,y) . 5.3.27 
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It is seen that a zero of P J is a fixed point of T^. We shall now show that 
(Py) ' (0.0) has a bounded inverse and the desired conclusions concerning the 
existence of a solution y may be deduced from the implicit function theorem. 

(For presentations of the implicit function theorem, see Kantorovich [K4] and 
Holtzman (HI]). The operator derivative is given as 

1 

(Py)' (0 ,0)v(t) = v(t) - J G J (t,s)(S + (|£)(0) - V(s) }v (s)ds 5.3.28 


or 


(Py)’ (0 ,0)v = [I-D J ]v 


5.3.29 


If the set J = {S + (|£-) (0) , M,N,} is boundary compatible, then [I-D J ] 1 is 

°y 

bounded and is given as (see Falb [FI]) 


[I-D J ]‘ 1 v = [I - R J ] v 


5.3.30 


where 


- r^J 

R J v = J G J (t,s){V(s) - S - (|£)(0)> v (s ) ds . 


5.3.31 


<v> gp 

All that now remains is to show that J = (S + (-r— ) (0) ,M,N) is boundary compatible. 

dy 

We have 


F(y) = 


<K*,u(p)) 

'(§7 ) 1 ( x > u (p))p . 


5.3.32 


and 


,8F, , .. 

^Hy) = 


cf|) (x,u(p)) 


- |^c||)'(x.»Cp)}p] 


(f|0 (x,u(p)) 


l^ulJj'Cx.uQpJJp] 


5.3.33 
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We have from (5.3.10) that (3tp/3x) (0,0) = 0 , and (3ip/3p) (x,u(p)) may be 
obtained as 


(|^( X , U (P)) = (|J)(|i)(x,u(p)) 


[(|f)(x,u(p)) - B] B ' 


which evaluated for y(-) = 0 yields 


(|»)( 0 , 0 ) . 0 . 


Defining the n x n matrix D(x,p) as 
D(x,p) = (|f) ' (x,u(p)) , 
it only remains to calculate 


|j[D(x,p)p] and [D(x,p)p] . 


If the n x n matrix Q(x,p) is defined as 


Q( x ,p) = [D(x,p)p] , 


Then the elements of the matrix are given as 

" 3D 

q k . (x,P) = ( 3F i) (x 'P^j 

1 j=l 1 


and then. 


Q(0,0) = 0 


Similarly if the n x n matrix T(x,p) = — [D(x,p)p] 

dp 


then the elements of the matrix are given as 


3D,. 


YkiCx.P) = X (sp^-Kx.PiPi + D(x,p) ki 

i-i 1 


5.3.34 


5.3.35 


5.3.36 


5.3.37 


5.3.38 


5.3.39 


5.3.41 


5.3.42 


80 



Since D(0,0) = 0 , then from (5.3.42) 


r (0,0) = o . 


5.3.43 


5.3.44 


5.3.45 

Then from the assumption that the set {A,B} is controllable and the result of 
Corollary 5.2.18, the set J = (S,M,N} is boundary compatible, and consequently 
the inverse is bounded. Hence for c in a neighborhood of the origin, a solution 
y exists to the TPBVP and the system is null controllable in a neighborhood of 
the origin as was to be proved. In addition, we note that the terminal state 
is not required to be the origin, but may be any point x^ in a neighborhood of 
the origin. 

The previous theorem does not specify the size of jt , the controllable 
region, only that the system is null controllable in a region near the origin. 

In addition, the condition that the linearized system be controllable about the 
origin is not necessary for nonlinear null controllability. The use of the 
contraction mapping theorem allows us to consider the domain of Xq,x^ such that 
a solution exists to the TPBVP, and moreover, the theorem is stated without 
specifying linearized controllability. As an example of the use of the contrac- 
tion mappings theorem for controllability investigation, the broad class of 
systems described as 


As a result. 


cff) (0) = 0 , 


and J is given simply as J = {S,M,N} where 


S = 


A -BB> 
0 -A' 
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x = f(x) + Bu 


5.3.46 


will be considered. 

Theorem 5.3.47. 

Consider the control process in R n 


x = f(x) + Bu in in R n x n 


5.3.48 


Let y q C * ) t> e an element of £>([0,1],R ) and let S = SCy^.r). Suppose that 
(i) J = {V,M,N} is a boundary compatible set, and (ii) there are real numbers 
n and a with n > 0 and 0 < a < 1 such that 


i) II T (y 0 )-y 0 



V 


. o . 


/eV, 


S) 


-BB’p 0 (s) + f(x 0 (s)) 







-x 0 (sn7 

x 0 (t) 


V (s) 

Po (s) J{ ds ' 

p 0 (t) - 



< n 5.3.49 


2) sup j 

||o6 * | 

j 

= sup_ 

sup 


y€ S 1 

ii y ii 

1 yes 

l|u|| 1 1 

J 

n 


G (t,s)f 


(||)Cx(s)) 


- ~ [c|^)’(x(s))p(s)] 


#'Cx(s)) 


3x LL ax-' 


-V(s)l u(s)ds 


5.3.50 


< a 


R ^r. 


5.3.51 


Then the CM sequence {y (’)} for the TPBVP based on y ^ and J converges uniformly 
to the unique solution y* in S and a control exists, u = -B’p*, to steer the 
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system (5.3.48) from x Q to the origin. 

Proof. Consider the optimization problem consisting of the system (5.3.48), the 
cost functional 

1 

J = j J <u(t), u (t )) dt, 5.3.52 

0 

and the boundary conditions 

x(0) = x Q , x(l) = 0 . 5.3.53 

Application of the Pontryagin principle to the posed optimization problem 
reduces the necessary conditions for optimality to the TPBVP 


' x 


- 

BB'p + f(x) 






p . 


- 

(M, 

'•ax - 1 

(x)p 






I 

0 


x(0) 

+ 

0 

0 

x(l)- 


x o 

0 

0 . 


p(0) 


.1 

0 

. PCD 


0 


5.3.54 


5.3.55 


For the boundary compatible set J = (V,M,N), the solution to the TPBVP may be 
written under certain smoothness conditions as 


x(t)' 

T J (y)(t) = H J (t) 

V 

p(t). 


.0 . 



- -BB'p(s) + f(x(s)) 

. -(|f)'(x(s))p(s) 


- V(s) 



5.3.56 
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Applying contraction mappings Theorem 3.4.14 to the operator (5.3.56) yields 
the conditions to be proved in Theorem 5.3.47. 

For the general system formulation, x = f(x,u), the canonical equations 
of 2.3.54 are used subject to the boundary conditions (5.3.55). Techniques for 
calculating the criteria of Theorem 5.3.47 will now be considered. 

5.4. Evaluation of Controllability Convergence Parameters 

In Section 4.6, the variables z^, z, and P(t) were defined such that coarse 
estimates for n and a were obtained as 

n = ilP(-)z n l| 5.4.1 

and 

a = flP(-)z II • 5.4.2 

In this section we shall consider the determination of z^, z, and P(t) for 
the controllability Theorem 5.3.47. In particular, the conditions for boundary 
compatibility and the use of simple V matrices and similarity transformations 
will be considered. 

The boundary value set for controllability problems is given as 



I 

0 ' 


0 

0 

M = 

0 

0 . 

N = 

.1 

0 


Assuming a general form for the 2n * 2n V matrix and the fundamental matrix 
® V (t,s), i.e.. 


■ V 11 

V 12 

v 

* (t,S) = 

« 11 (t,s) 


■ v 21 

V 22 


fl 21 (t,s) 

« 22 (t,s) . 


The matrix [M+N4> V (1 ,0)] is obtained as 
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[M+N$ (1,0)] = 


0 


5.4.5 


l Q n (l,0) fl 12 (l,0) J . 

This yields det[M+N$ V (l ,0)] = det C 1 >0) 3 • Hence the condition for boundary 
compatibility reduces to the nonsingularity of fl^Cl,®)* In passing, it is seen 
that neither the zero matrix nor a diagonal matrix (nor a modified diagonal) may 
be used in the integral representation. If det (1 ,0) ] ^ 0 , the inverse of 
[M+N0 V (1 ,0) ] is given as 


[M+N4> V (1 ,0) ] -1 = 


I 

-fi^(l,0) fl n (l,0) 


0 

fl'^Cl.O) . 


5.4.6 


and then the core matrices of the Green's function are given as 


[M+N0 V (1,0)] _1 M 


-njj (i,o) n n (i,o) o 


and 


[M+N* V (1,0.)] _1 N* V (1,0) 


0 1 


fi 12 (1 ’ 0) n il C1 ' 05 1 


5.4.7 


5.4.8 


Since the boundary value set (5.4.3) disallows the use of particularly simple 
V matrices, we shall consider an approximate technique for calculating G J (t,s) 
and P(t) for V matrices of general structure. 

Example 5.4.9. 

Using the canonical transformation 


D = A _1 VA , 


5.4.10 
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the Green's function matrices are given as 


Gj(t,s) = {A$ D (t,0)A' 1 }{[M+N<I> V (l,0)]' 1 M}{A'5> D (0 ) s)A' 1 } 5 

Gjj (t ,s) = -{A$ D (t ) 0)r 1 H[M+N4 V (l l O)]'V(l,0)HA4 D (0,s)A' 1 } S 

for the boundary compatible set J = (V,M,N). From (5.4.7) and (5.4.8), the 
center bracketed terms in (5.4.11), (5.4.12) are given as 

[M+N$ V (1,0)]' 1 M = 

and 



5. 
5. 

Assuming that 4>°(t,s) represents the primary magnitude characteristics of 
A4> D (t,s)A" 1 , approximations for Gj(t,s) and Gjj(t,s) are. formed as 

G ^ (t,s)2 (t ,0) [M+N$ V (1,0)] " 1 M4 >D (0 ,s) 5. 

and 

GjjCtjS) '■'1 -$°(t,0) [M+N<S> V (l,0)]”^N4 V (l,0)5> D (0,s) . 5. 


where 


[M+N$ V (1,0)]‘ 1 N4 V (1,0) = 


A = Ojjd.O) 


0 0 ‘ 

-A I, 


Following the discussion in Section 4.6, V is chosen such that $ D (t,s) is 
diagonal or modified diagonal and may be represented as 


$ D (t,s) = 


K 11 


(t,s) 0 


0 ^22 (t > s) 



.4.11 

.4.12 

4.13 

4.14 

4.15 

4.16 

4.17 

4.18 
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Using this form of 4 (t,s) in (5.4.16) and (5.4.17), the approximations for 
Gj(t,s) and Gjj(t,s) are given as 


Gj (t,s) 


'11 


(t,s) 


and 


(t*s) ~ 


^(t.O) (0» s ) 0 


i° 2 (t,0) A^O.s) ^(t.s) 


5.4.19 


5.4.20 


The approximation for P(t) is then given as 


P(t) 


L. 

/' 


11 


(t.s) ds 


J. 

/ 

L 0 


l*2 2 (t«°) Mn(0,s)|ds 


'll’ 


x 

f k22 (t * 


s) | ds 


5.4.21 


where 


,-l 


A = -n 12 (i,o) 0 n (l,0) 


5.4.22 


For D a diagonal matrix, P(t) given by (5.4.21) becomes a particularly simple 
form. Several variables are now defined which will be used with P(t) to calculate 
estimates for n and a in Theorem 5.3.47. 

Definition 5.4.22. 

Using the boundary compatible initial estimate, y^(t) = H^(t)c, define the 
2n vector z Q as 
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|-BB'p 0 (s) - V n x 0 (s) - V 12 p 0 (s) + f(x 0 (s)) 


z = sup 

0 s I _V 21 X 0 ts) - V 99 Pn( s ) - (3f/8x)' (x n (s))p n (s) | 


22^0 


'O' ■’"‘'O' 


5.4.23 


A conservative estimate for n is given as 


= l|P(Oz f 


5.4.24 


Definition 5.4.25. 


Define the real numbers v. .. £. .. and a. . as 

ij iJ ij 


Vi j = Th ll ( 33T )(x) ' V ll..l \> i’i* 1 ’* 

3 ij 


5.4.26 


and 


= I C- BB ’)ii - V 12 I ; ij=l.n 

J ij 


I ( ^% Kx) Pkl I • 


5.4.27 


5.4.28 


and define the n vectors Zj. and z to be composed of the elements 


r -r ~ ) ' (V. . + £. .) 

Ij A—rf IJ XJ 


j=l 


5.4.29 


and 


hi. = £ ( a ij + V 
1 J=i 


5.4.30 


Then the 2n vector z defined as 

Z I 
Z II 


5.4.31 
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may be used with P(t) to obtain a coarse estimate for a as 

a = ||P(Oz! • 5.4.32 

Numerical evaluation of the convergence criteria is presented for various 
examples in Chapter 6. 
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CHAPTER 6 


iNUMEKlCAL EXAMPLES 


6.1. Introduction 

We examine the regulation and control of several nonlinear systems to 
demonstrate the usefulness of contraction mappings and to illustrate the practical 
applications of the major theorems. There are many well known and very powerful 
iterative methods for the solution of optimal control problems. However, practical 
convergence criteria are few and far between. In this chapter we demonstrate 
that general results may be obtained via the application of the contraction 
mappings convergence theorem. In addition, the practical application of the 
contraction mapping algorithm demonstrates that in many cases it is an efficient, 
straightforward technique for the solution of optimal control problems. The 
examples demonstrate that practical application has a much broader range than 
the theoretical results might imply. This is primarily due to the coarse 
estimates which are used to evaluate the convergence parameters. 

The first example to be considered is the regulation of the well known 
Van der Pol equation. As an illustrative exercise, both contraction mappings 
and modified contraction mappings are applied to this problem. The results 
obtained are compared with previously published data. The second example 
begins a two part sequence investigating the null controllability of nonlinear 
systems. The first member of the sequence is a particularly simple system 
which serves to introduce bounded control problems. The final example of the 
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chapter considers the null controllability of the pitching motion of a satellite 
with bounded control thrust. 


6.2. Van der Pol's Equation 

In this rather long example, we consider in detail many of the concepts 
essential to the contraction mappings theory. In particular, the choice of 
the boundary compatible set J and the calculation of the convergence parameters 
will be investigated closely. The system to be considered is the driven, second 
order nonlinear oscillator studied by Van der Pol . 


-x + £.(l-x 1 )x 2 + u . 


6 . 2.1 


The cost functional to be minimized is taken from Bullock [B6] as 
5 

J = J f [Xj(t) + x 2 (t) + U 2 (t)]dt , 6.2.2 

0 

and the boundary conditions for the optimal regulator problem are given as 


x^O) = 1.0 x 1 (5) unspecified 

x 2 (0) = 0.0 x 2 (5) unspecified . 6.2.3 


From Example 2.3.10, the necessary conditions 


for optimality reduce to the TPBVP 


y = Sy + f(y) 6.2.4 

Ky(0) + Ly(l) = c 

where 


' 0 s o o' 


0 



2. 

-5 0 0 -5 


5(l-x )x 


f(y) = £ 

1 £ 

-5 0 0 5 


10x l x 2P 2 

. 0 -5 -5 0. 


-5(1- x 2 )p 2 
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and 



'l 

0 

0 

o' 


' 0 

0 

0 

o' 


1' 


0 

1 

0 

0 


0 

0 

0 

0 


0 

K = 

0 

0 

0 

0 

L = 

0 

0 

1 

0 

c = 

0 


.0 

0 

0 

0. 


. 0 

0 

0 

1. 


0. 


Using a boundary compatible set J = (V(t), M, N), (6.2.4) may be expressed in 
integral form as 


y (t) = H J (t){c-Ky(0) - Ly Cl) + My.(O) + Ny Cl) ) 

1 

♦ yV(t,s){Sy(s) + f(yCs)) - V(s)y(s)}ds 
0 


6.2.7 


The iterative solution of (6.2.7) by contraction mappings is now considered. 
We begin with the selection of the boundary compatible set J = (V(t),M,N). 

Since the boundary conditions of (6.2.4) are linear, the natural choice 
for the matrices M and N are M=K, N=L. If the initial estimate of the solution 
is then chosen as H^Ctjc, every member of the contraction mapping sequence 
satisfies the boundary conditions. As indicated previously, it is often 
advantageous to choose the matrix V in such a way that (S + (3f/3y) (y) -V(s) } 
is small. Following this guideline generally requires inclusion of time 
varying functions in the V matrix, thus complicating the convergence analysis. 
However, if £ is small in (6.2.5), an acceptable choice for V is simply the 
linear part -of (6.2.4), i.e., V = S. For larger values of £ , it may become 
necessary to include an effect of the nonlinearity in the choice of the V matrix, 
but for now, consider V to be chosen as 
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0 5 0 0 


-5 

-5 


0 

0 


L 0 -5 


0 -5 

0 5 

-5 0 


6 . 2.8 


The variables P(t), z Q , and z, defined in Chapter 4 for use in convergence 
analysis, will now be obtained for this example. From (4.2.9), (4.2.11), (6.2.5) 
and (6.2.8), the vectors z Q and z are 


z 0 = sup 


1 5 (1-Xj(t) 0 )x 2 (t) 0 | 

1 10 x 1 (t) 0 x 2 (t) 0 p 2 (t) 0 | 


|S (l-xj(t) 0 )p 2 (t) 0 | 


6.2.9 


and 


r 

0 



1 2 x x (t)x 2 (t)| + 

Id-xj(t))| f 

£< 

z = sup sup_ \ 

1 2 x 2 (t)p 2 (t)| + 

1 2 x 1 (t)p 2 (t)| + 1 2 Xl (t)x 2 (t)| / 

t y€ S 


1 

( 

1 2 X j (t)p 2 (t) | + 

(l-x^(t))| JJ 


for V given by (6.2.8), the calculation and structure of the fundamental matrix 
is somewhat involved, thus complicating the calculation of P(t). Since the 
characteristic roots of V are two pairs of complex conjugates, the techniques of 
Example 4.4.24 are useful for evaluating the P(t) matrix. The canonical 
transformation 

D = A _1 VA 6.2.11 
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transforms the V matrix into the block diagonal form 


0 0 

Oj 0 0 . 

0 0 a 2 “2 

0 0 -a> 2 a 2 


6 . 2.12 


where =-3.35, = 4.91, and a 2 = 3.35, to 2 = 4.91. It is then straightforward 

to calculate the matrix P(t) from (4.7.1) and (4.7.2). 

The parameter £ is included in the example so that the general case may 
be considered. In particular we are interested in determining the range of £ 
for which the contraction mappings theorem is valid. Before proceeding with 
the convergence analysis, the sphere S(yg,r) must be defined. The initial 
estimate of the solution, y.(*) is taken to be the boundary compatible initial 
estimate, i.e,, the solution to the linear TPBVP 


y = Vy My(0) + Ny(l) = c 

or 6.2.13 

y (t) = H J (t) c. 


(It should be noticed that this choice for y^ does not require additional 
computation since the terms are necessary for the CM algorithm.) With this 
choice of y , the radius of S is taken as r = 0.1. This sphere S(y 0 ,r) is 
illustrated in Figure 6.1. 
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Figure 6.1 The Sphere S(y 0 ,r) 

From (6.2.9) and (6.2.10), the vectors z Q and z are calculated as 



Conservative estimates for the convergence parameters r\ and o are obtained as 

n = sup(P(t)z 0 ) = 2.1 £ 6.2.15 

t 

and 


a = sup{P(t)z) = 6.7£ 


6.2.16 



Using (6.2.15) and (6.2.16), the requirements of Theorem 3,4.14 are specified as 


and 


6.7 £ < 1 


.1 £ 

1-6. 7£ 


0.1 


6.2.17 

6.2.18 


Analysis of (6.2.17) and (6.2.18) shows that for £ <_ 0.034 the convergence 
conditions of the theorem are satisfied. 

The case for £ = 1.0 is treated in the paper "A Second-Order Feedback 
Method for Optimal Control Computations", by Bullock and Franklin [B6] . In 
the paper, the optimization problem presented by (6. 1,2, 3) is solved by the 
techniques of steepest descent and second variation. We now consider the 
application of contraction mappings to (6.2.4) with £ = 1.0. Again the V matrix 
is chosen V = S. Now rather than taking y^ as the solution to the linear TPBVP, 
yg(t) shall be given by the fifth iteration of the CM algorithm begun with the 
initial guess If 1 (t)c. This choice for y Q is made so that the region S(y Q ,r) is 
more likely to include the solution y(t) to the nonlinear TPBVP. We again take 
r = 0.1. This sphere is illustrated in Figure 6.2. 

We shall first determine the convergence rate factor a. Taking the supremum 
over S, z is given as 


0 . 0 " 

6.1 

12.9 

. 8.1J , 
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Figure 6.2 The Sphere S(y Q> r) 


and a coarse estimate for a is 

a = sup{P(t)z) = 6.8 6.2.19 

t 

With a > 1, the conditions of the theorem are not satisfied and convergence is 
not guaranteed by the theorem. However, these theoretical results are only 
guidelines for the practical application of contraction mappings. In fact. 
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the CM algorithm reduced the convergence norm (given as sup sup|y^ j(t)-y^ n (t)|) 

ie p t ’ ’ 

-4 

to 10 in fourteen iterations. In order to compare these results with those 
presented in [B6] , we note that the norm of the cost function (given as |j J n l) 

was reduced to 10 ^ in fourteen iterations by the CM algorithm. In [B6], the 
computed cost agreed with the optimal in only two significant figures after 
eighteen iterations for the steepest descent procedure. The more complicated 
second order technique obtained five place accuracy in the cost after five 
iterations . 

Using the results of Section 3.5, we now consider a technique which is 

often effective in reducing a and the number of iterations required by the CM 

algorithm. In this approach, a more complicated boundary compatible set 

J = (W(t), M, N} is used in the integral representation. The matrix W(t) is 

designed to include time varying terms attempting to model the effects of the 

2 2 

nonlinearity. For example, model [1 - Xj(t)] as [1 - (1-t) ] and select the 
W(t ) matrix as 


W(t) 


0 5 0 0 

-5 5£[l-(l-t) 2 ] 0 -5 


-5 0 0 5 

0 -5 -5 -5 £ [ 1- (1-t) 2 ] 


6 . 2.20 


However, using the equivalence relation (3.5.29), we have 


dfy) ■ U-uirVw -.>4 y] 


6 . 2.21 


for the boundary compatible sets J = (V,M,N) and J = {W(t),M,N}. Hence the Green's 
function may be calculated using the simpler set J = {V,M,N} where V is given by 
(6.2.8) and P(t) is calculated using the D matrix (6.2.12). We previously found 
that with V given by (6.2.8), the conditions of the contraction mappings theorem 
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are satisfied for d < 0.034. A similar analysis is now done for J = (W(t),M,N). 


The vectors z Q and z axe given as 


z Q = sup 


|5x 2 (t) 0 [l-Xj(t) 0 ) - (1-Q-t) 2 )] 


1 10Xj (t) 0 x 2 Ct) 0 p 2 Ct) 0 l 
[|5p 2 Ct) 0 [l-x 2 (t) 0 ) - (l-(l-t) 2 )]|J 


6 . 2.22 


and 


z = sup 
yeS 


|2x 1 (t)x 2 (t)| ♦ | (l-x 2 (t)) - (l-(l-t) 2 )| 
|2x 2 (t)p 2 (t)| + |2x 1 (t)p 2 (t) | + |2 Xl (t)x 2 (t) | 
|2x.(t)p,(t)| ♦ |(l-x?(t)) - (l-(l-t) 2 )| 


6.2.23 


Now using J = (W(t), M, N} with y Q (t) = H J (t)c, r = 0.1, we find following 
Example 3.5.25 and (6.2.22), (6.2.23) that conservative values for the convergence 
parameters are 


and 


n = sup {P(t)z Q } = 0.64£ 


a = sup (P(t)z) = 5. Of 
t 


The requirements of Theorem 3.4.14 are then 


and 


5. Of < 1 


0.646 „ 

1-5.06 - 


0 . 1 . 


6.2.24 

6.2.25 
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Analysis of (6.2.24), (6.2.25) shows that the convergence conditions of the 
theorem are satisfied for 

i < 0.092 6.2.26 

a three fold increase over the previous value. These results are guidelines, 
but indicate the improvement due to use of the better designed, though more 
complicated, W(t) matrix. 

Using the boundary compatible set 3” = {W(t),M,N} for & = 1.0, a conservative 
value for a is a = 5.1, an improvement over (6.2.14), but again violating the 
theoretical specifications. However, the practical application of the CM 
algorithm reduced the convergence norm to 10 ^ in eight iterations, a significant 
improvement over the algorithm using J = (V,M,N). The iterative sequence for 
the control function is shown in Figure 6.3. 



Figure 6.3 Control Iterations 
(Numbers indicate iteration sequence.) 


101 



A comparison of the convergence behavior for J and J is shown in Figure 6.4. 

NORM 



Figure 6.4 Comparison of Performance for Contraction 

Mappings and Modified Contraction Mapptions . 

6.3. Null Controllability with Bounded Control 

The first example of system null controllability involves a simple linear 
system with bounded input control. The example is included primarily as an 
introduction to the techniques of dealing with a bounded control. Consider 
the system 

x = Ax + Bu 6.3.1 
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where 



0 

f 


'0 

A = 



B = 



.0 

0 . 


.1 


and the control magnitude is constrained to satisfy 

j u (t ) [ <^ 1 , 0 <_ t <_ T . 6.3.3 

The initial conditions are 

XjCO) = 1 , x 2 (0) = 1 , 

and the final state of the system is required to be the origin, i.e., 

Xj(t) = 0 , x 2 (T) = 0 , 6.3.5 

where T is a prescribed fixed terminal time. 

The linear system (6.3.1) is clearly controllable since rank [B,AB] = 2. 

However there do exist combinations of T and x Q such that the system cannot 

be driven to the origin by the bounded control in time T. We shall investigate 

the null controllability of this system by considering the optimization problem 

composed of the system (6.4.1), the cost functional 
T 

J=y /u 2 (t)dt, 6.3.6 

0 

and the boundary conditions (6.3.4) , (6.3.5) . 

Analytical investigation of this optimization problem yields the 
information that the minimum time required fox the system to be driven from 
(1,1) to the origin is 1 and at this value, the H-minimal control is 

bang-bang. As T is increased from 1 + VT, the optimal control becomes a 
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saturating function, and wlien T is sufficiently great, the H-minimal control 
never saturates, i.e., it never takes on its maximum allowable magnitude. These 
points concerning null controllability are now illustrated by applying contraction 
mappings to the TPBVP associated with the posed optimization problem. 

Application of the minimum principle and a change of time variable transforms 
the optimization problem into the TPBVP 



where SAT(-) is defined in (2.2.17), and where the time variable has been 
changed so that t = as where s 6 [0,1] and a = T. [(•) now indicates differentia- 
tion with respect to s.] We shall consider the case with a = 5.0. 

Using the boundary compatible set J = (V,M,N), the solution to (6.3.9), 

if it exists, may be written as 

1 

y (t) = H J (t)c ♦ / G J (t,s){f(y(s)) - V(s) y (s) }ds . 6.3.10 

0 
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From Corollary 5.2.18, the V matrix given by 


V = 


A 

0 


-BB' 
-A' . 


6.3.11 


is boundary compatible with M and N given in (6.3.8). Using V specified in 
(6.3.11), (6.3.10) may be written explicitly as 


y(t)=T J (y) (t)=H J (t)c + 


JL 

f 


G J (t,s) 


0 

ap 2 (s) -aSAT{p 2 (s) } 
0 
0 


ds 


6.3.12 


We shall now investigate the convergence conditions for the contraction mappings 
algorithm when applied to this non-differentiable TPBVP. Instead of deriving 
conditions satisfied by the Frechet derivative, we shall be concerned rather with 
conditions on the Lipschitz norm of the operator T*. The initial estimate of 
the solution and the center of the region S(y Q ,r) is taken to be H J (t)c. To 
complete the definition of S, the radius is set as r=0.2. This region is illus- 
trated in Figure 6.5. 

Values for || T J (y Q ) -y Q || and the Lipschitz norm 0t J 0 - must now be 
calculated. From (6.3.12), it is seen that the nonlinearity is contained in 

only the second component of the forcing function. Hence we shall investigate 

the Lipschitz norm of the operators 

1 

T^(u) = yv 2 (t,s)[au(s) - aSAT{u(s) }]ds 6.3.13 

0 

where G^is the i,2 element of the Green's matrix G^(t,s). The Lipschitz 
norm is formally defined as 
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Figure 6.5 The Sphere S(y Q ,r) 


D T i D § = sup - 

u,v 6 S 
u* v 


H T i(u) - T.(V)|| 
llu - v|| 


6.3.14 


and for the operator in (6; 4. 13) 


|| T. (u)-T.(v)|| l Jo < 2 (t,s) [a(u (s)-vfs) ) + aSAT{v(s) }-aSAT{u(s) }]ds j| 


I|u-V|| 


IlufO - v(.)ll 


6.3.15 

or 
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1 

j|T i fu)-T i Cv)|| a sup J | 2 tt: , s ) { [u(s) -v(s)}+[SAT{v (s) }-SAT{u(s) } j } | ds 


< t 0 


u - v 


sup |u(p)-v(p) 


6.3.16 


Now noting 


sup|u(p)-v(p) | > |u(s)-v(s)| 


6.3.17 


and 


sup |u(p) -v(p) | >_ | SAT{u(s) }-SAT {v(s)}| , 


6.3.18 


we have 


|| T i Cu)- T i Cv)|J 

I|U - V|| 


< a sup 
t 


/|< 2 (t,S) 


0 


u(s)-v(s) SAT{v(s) }-SAT{u(s) } )i , 
|u(s)-v(.s)| + |u(s)-v(s) | jl dS 

6.3.19 


It may be shown that with r = 0.2 


( u(s)-v(s) + SAT{v (s) }-SAT{u(s) } 

||u(s)-v(s)| |u(s)-v(s)[ 


The required Lipschitz norm may now be evaluated in several ways, each of varying 
degrees of accuracy. We now consider one of the more accurate techniques. 

Note that as a result of the choice of S(yg,r), saturation can only occur 
for 0 < s < 0.25. We may now write 


HV u) -V v) li 

||u - v|| 



6.3.21 


which may be approximated as 
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6.3.22 


Ql^D = < sup sup {2a(.25)P (t)} = S/16 < 1 

b - i t 12 


where 


Now 


P i2 (t) = / I G 1 Ct,s)| d s + f \G J U (t,s)|ds 

0 1,2 t 1,2 

determining ||t J (y Q ) - y Q |l , we have 


(y 0 ) - y 0i = f G^ 2 (t,s) [ap 2 (s) 0 - aSAT(p 2 (s) 0 }]ds 


From y (•) 
7 o 


sup( | ap 2 (t) Q - aSAT{p 2 (t) Q } [ } <_ 0.04a , 


and then 


.25 


|T J (y 0 )-y 0 ll <_ sup sup{0.04a f jG. 2 (t,s) jds} 
"’it o' 

< (0.04a) (0.25)sup sup (P. 2 (t)} = . 


i t 


Taking conservative values a = S/16, n = 1/60 , 


1-a 


0.1 < r = 0.2 


so that the theoretical application of contraction mappings is successful. 
Hence a solution exists to the TPBVP and a control exists to accomplish the 
desired transfer. These concepts will now be applied to a nonlinear system. 


3.23 

3.24 

3. 25 

3.26 

3.27 
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6. 4. Controllability of Satellite Pitch Motion 

The pitch motion of a satellite in circular orbit can be described by the 
normalized differential equation 


X- 


f X- 

1 

- 

2 



-sin x, +u 

2 J 


1 


whenever a principal axis of the satellite remains normal to the orbit plane 
[Rl] . The controlling torque u(t) is bounded (|u(t)| <_ 1) and x^(t) is twice 
the pitch coordinate. In investigating the null controllability of the system, 
our goal is to find an acceptable control u(t) which zeros the pitch and pitch 
rate in a prescribed fixed time T . 

In a neighborhood of the origin, system (6.4.1) behaves as 


x = Ax + Bu 6.4.2 

where 



' 0 1 


0' 

A = (8 f/3x) (0,0) = 


B = (3f/3u) (0,0) = 



.-l o. 


1 . 


For the linear system (6.4.2), rank [B,AB] = 2, and from Theorem 5.3.3, the 
nonlinear system (6.4.1) is controllable in a region of the origin. As in the 
previous example, null controllability is investigated by considering the optimiza- 
tion problem consisting of the system (6.4.1), the specified boundary conditions, 
and the cost functional 


T 

J = \ f u 2 (t)dt 
0 


6.4.4 


109 



Application of the minimum principle and a change of time variable trans- 
forms the optimization problem into a TPBVP of the form 


Y = Sy + f(y) 
My(0) + Ny(l) = c 

where 



0 

a 

0 

0 



r 


0 

1 


0 

0 

0 

0 







s = 

0 

0 

0 

0 

f(y) = 

-a 

sin x^-aSAT{p 2 ) 


.0 

0 

-a 

0. 




ap 2 cosx^ 









- 


0 

- 


1 

0 

0 

0 


0 

0 

0 

0 ‘ 


o 
!— < 

X 

M = 

0 

1 

0 

0 

N = 

0 

0 

0 

0 

c = 

% 


0 

0 

0 

0 


1 

0 

0 

0 


0 


.0 

0 

0 

0 


.0 

1 

0 

0. 


0 


6.4.5 


6.4.6 


6.4.7 


where the differentiation is now with respect to s where t = as , s£ [0,1]. 
If a solution to (6.4.5) exists, it may be represented as 


1 

y(t) = H J (t)c + yG J (t,s){Sy(s) + f (y(s)) - V(s)y(s)}ds 6.4.8 

0 


where J = (V(t), M, N} is a boundary compatible set. Since the linear system 
(6.4.2) is controllable. Corollary 5.2.18 states that the 2n x 2n V matrix 
given as 


V = 


A 

io 


-BB' 
-A' . 


is boundary compatible with M and N given by (6.4.7). 
yields 


6.4.9 

Choosing V from (6.4.9) 
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V = 


0 a 0 <5 

-a 0 0 -a 

0 0 0 a 

.0 0 -a 0. 


6.4.10 


We shall now investigate the null controllability of the system for various 
initial conditions and time intervals. 


Example 6.4.11. 

Consider the initial condition for the system to be 60° for the pitch angle 
and zero for the pitch rate. It is desired to regulate the system to zero in 
one half period, i.e., T = it. The initial estimate of the solution and the 
center of the sphere S is taken to be H^(t)c. The region SCy^.r) is defined by 
setting r = 0.1 and is illustrated in Figure 6.6. 

It is seen that for this S(yg, r) that |p 2 | is always less than one so that 
saturation never occurs. Hence the forcing function for (6.4.8) may be con- 
sidered as 

r o i 


F(y) = Sy + f(y) - Vy = 


-a(sin Xj - xp 
a(cosx 1 - l)p 2 
0 


Estimates for the convergence parameters are calculated using the variables 
defined as 


1 

P(t) = J |G J (t,s)|ds 6.4.13 

0 


z 0 = SU P { I F i. (y 0 Cs) D | > 6.4.14 

i s 
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Figure 6.6 


The Sphere S(y QJ r) for T = ir. 


and 

( P l/ 3F '\ 

Vi"Y £ 1(557) o-csj) 

y£b < j=i 3 

where F(y) is given by (6.4.12). The matrix V given by (6.4.10) may be transformed 
into the canonical matrix D given as 



2 


CJj 'Mj 0 0 

0 ^ 0 0 


0 0 

0 0 -u 


°2 w 2 


2 °2 


6.4.16 


where cr^ = 0, io^ = a, = 0.05, w 2 = a. From (6.4.14), the matrix P(t) may then 
be obtained by integration of the expression 


P(t) = y*|{A4 D (t,0)A' 1 H[M+N4 V (l,0)]' 1 M}{A4 D (0,s)A' 1 }|ds 


y* |{A$ D (t,0)A' 1 H[M+N0 V (l,0)]' 1 N4 V (l,0)}{A$ D (0,s)A" 1 }|ds 


Using (6.4.12), (6.4.15), and taking the supremum over S yields 

0 

| a(cosx^-l) | 

|ap 2 sin Xj| + |a(cos x 1 -l) | 

0 


z = sup_ 
y6 S 


1 

0.0 • 


0.418 


0.653 


.0.0 


6.4.17 


6.4.18 


The vector z Q is found from (6.4.12) and (6.4.14) as 


0.0 

0.095 

0.028 

0.0 


6.4.19 
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Using (6.4.17), (6.4.18), and (6.4.19), conservative estimates for the convergence 
parameters n and a are found as 

n = sup (P(t)z } = 0.047 6.4.20 

t U 

a = sup (P(t)z) = 0.38 6.4.21 

Now testing n/(l-a) < r , we have 


J2- = 0.075 < r = 0.1 

1-a 


6.4.22 


Hence the convergence conditions of the contraction mappings theorem are 
satisfied and the theoretical application of contraction mappings is successful. 
Moreover, a solution exists to the TPBVP and a control exists to accomplish the 
desired transfer. 

Example 6.4.23 

Consider the initial condition for the system to be 60° for the pitch angle 
and zero for the pitch rate. It is desired to regulate the system to zero in one 
quarter period, i.e., T = m/2. The initial estimate of the solution and the 
center of the sphere S is taken to be If 1 (t)c. The region S(y 0 ,r) is defined by 
setting r = 0.1. This region is illustrated in Figure 6.7. 

It is seen that S(y Q ,r) contains a saturating region for p 2> Hence the 
forcing function for (6.4.8) must be considered as 


F(y) = Sy + f (y) - Vy 


0 

-a(sin x^xp - a(SAT{p 2 )-p 2 ) 
ap 2 (cosxj-l) 

0 


6.4.24 
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Figure 6.7 The Sphere S(y Q ,r) for T = u/2. 


As in Section 6.3, the Lipschitz norm of the operator (y ) , not the Frechet 
derivative, must be investigated. For SCy^.r), the Lipschitz condition for 
F(y) is given as 


| F (y) -F (y ' ) | < 


' 0 

0 

0 

0 


iv x 

il- 

0.39 

0 

0 

3.14 


X 

i 

CM 

X 

2 I 

0.75 

0 

0 

0.34 


Ipj-p 

il 

.0 

0 

0 

0 


Ip 2 -p 

2 1 ■ 


6.4.25 
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or 

||F(y) - F (y ' ) || 1 C |]y - y'|| . 


6.4.26 


Using P(t) from (6.4.17) and defining the 2n vector z as composed of the elements 



6.4.27 


a conservative estimate for a is (|P( , ) Z II • However, because the saturation 

occurs only over a short interval of S(yg, r) , this estimate would tend to be 
quite inaccurate. Hence we deal with the saturation effect separately. Now let 


0 


fo 


z 

1 


0.34 

1.14 


and z^ 


3.14 

0 


0 


10 


6.4.28 


where z^ arises from the differentiable part and z 2 from the saturating effect. 
Then as in Section 6.3 , 


a = sup (P(t)z } + sup 
t t 


, 1 '° 

f |G J (t,s)z 2 |ds 
( 0.. 85 


6.4.29 


or 

a = sup {P(t)z 1 > + sup{(3.14) (0.15)P i2 (t)} . 


6.4.30 


Using the values of P(t) from (6.4.17), a is evaluated as 


a = 1.56 


6.4.31 
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Hence the requirement that a < 1.0 is violated, and convergence is not guaranteed 
However, as indicated previously, these coarse estimates are used as guidelines 
for the practical application of contraction mappings. Indeed, the CM algorithm 
reduced the convergence norm to 10 ^ in ten iterations. Figure 6.8 illustrates 
the state and control history. 



6.8 State and Control History 
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CHAPTER 7 


PRELIMINARY STUDY ON THE DYNAMICS OF 
DRUG USAGE WITHIN A COMMUNITY 


7.1. Introduction 

The modeling of complex socio-economic systems has recently received con- 
siderable attention. Arising initially as an aid to management decision making, 
[F7] , [R2], system modeling is now applied to many systems of public concern 
[F 5 ] , [R3] . The primary objective of the modeling effort is the formulation of 
improved administrative control policies. Typically, once a model is developed, 
the process of designing improved policies is largely a trial and error process. 
That is, the behavior of the system is first simulated with the model using one 
control policy and then another. The simulation results are then compared to 
determine which policy yielded the "best" behavior, clearly an inexact and 
inefficient technique of analysis. 

In this chapter we consider the feasibility of applying the systematic 
techniques of optimal control theory to the determination of policies for 
social systems. Specifically, a dynamic model attempting to represent the 
causal, feedback structure of community drug usage is developed. Then using 
optimization theory, we attempt to gain insight into how a community might best 
respond to a rapidly growing heroin addiction problem. The initial phase of 
the study is the creation of a dynamic model which reflects the modes of behavior 
of the system being investigated. 
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7.2. Development of a Dynamic Model 

The development of a model for a complex system such as drug usage is in 
itself a major effort. The subtle interrelationships and multi-feedback loops 
are often difficult to conceptualize. Likewise, the determination of various 
parameters within the model is a difficult task involving much data analysis. 

The main thrust of this chapter is not in the modeling direction. Rather, we 
shall develop a simple model which hopefully reflects inpart the basic behavior 
of a very complex system. Similarly, parameter values are chosen after discussions 
and readings and are believed to be reasonable. In this spirit, the development 
of the model is begun. (For the development of a more comprehensive model, 
see Roberts [R3]). 

The model concerns itself with three groups of people within the community. 
These groups represent the three levels of drug usage which will be considered in 
the model. These three pools of people are: 

i) potential drug users 

ii) drug users 

iii) heroin addicts 

Of course, much finer lines may be drawn, but these three are sufficient for this 
study. The dynamical nature of this problem is reflected in the constantly 
changing population of each level and the inherent relationships between these 
changes. The multiple interrelationships are often difficult to conceptualize, 
but are critical to the feedback, multiloop structure of the system. Figure 7.1 
represents how one might initially conceive this system as simply involving 
transitions of people from various stages of drug usage. 

In Figure 7.1, the double lines represent the flow of people between levels and 
the values controlling these flows are determined by the variables alongside. 
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Note that. in reality there are feedback paths from the drug users level and 
the heroin addicts level back to the community response. However, we shall 
be attempting to determine how a community might best respond rather than 
modeling the present reaction. 

Let us first consider the flow between potential drug users and drug users. 

In this study, "Potential Users" will represent the community population between 
the ages of ten and thirty who are not using drugs illegally. The next level of 
drug usage, "Drug Users," represents that group of people who occasionally 
participate in the illegal use of drugs, but who are not addicted to heroin. 

The flow between these levels is determined by the drug education program, the 
police effort, the number of potential users, and the number of drug users. Of 
these four variables, the number of drug users might be considered the dominant. 
This is simply due to the fact that the users tend to share their supply, turn-on 
their friends, and in general, tend to increase their numbers. The level of the 
drug education program and the fear of arrest may tend to deter some potential 
users, but these are not the dominant effects. The flow rate from potential users 
to users depends on the number of potential users in the sense of availability, 
i.e., if there are few potential users remaining, the inflow into drug users will 
wither, and, conversely, if there are many potential users, the self induced 
growth rate of drug usage is unimpeded. Some drug users revert back to potential 
users through the efforts of police and education, but this is considered a minor 
effect . 

The flow from drug users to addicts is of the same form as the flow from 
potential users to users. Again, education and police effort tend to deter the 
flow and a self growth rate is again present via the number of addicts. The flow 
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from Drug Users to Heroin Addicts simply reflects the fact that most addicts 
previously used "soft" drugs; it is not a causal indication. Addicts are removed 
from the street primarily by police action which is a result of the community 
response to the number of addicts manifested by the rising crime rate. The drug 
education program and police effort are created by community spending for these 
programs. In this model, the community spending for police and education are 
considered as the two control variables to be determined. 

For simulation and optimization studies, the general description of the model 
must be transformed into a system of equations characterizing the dynamics of the 
system. A convenient procedure for developing equations describing the dynamics 
of a general system is the DYNAMO format [P3]. Developed by the Industrial 
Dynamics Group at the Sloan School, M.I.T., DYNAMO is both a simulation language 
and a discrete equation representation for the system dynamics. We now develop 
the DYNAMO equations which describe the dynamics of the drug usage model. 

As indicated in the general description of the system, the number of drug 
users determines the nominal growth rate of drug usage, i.e,, the "recruitment" 
rate. This is represented as 

NGRU.K = DU.K 7.2.1 

where 

NGRU Nominal Growth Rate of Drug Usage 

.K a postscript indicating that NGRU.K 
refers to nominal growth rate at 
the present time K 
DU Drug Users 

AODC a constant determining the growth rate. 
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The availability of potential users is included as a multiplier of the nominal 
growth rate and is a function of the difference between the initial number of 
potential users and the present number of drug users. The nonlinear relationship 
has the general form illustrated in Figure 7.3 where APUM is the availability of 
potential users multiplier and IPU is the initial number of potential users. 


APUM 



Figure 7.3 Availability of Potential Users Multiplier 


The total flow from potential users to drug users is then given as 

GRU.KL = (APUM.K) (NGRU.K) 7.2.2 

where 
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GRU Growth Rate of Usage 

. KL postscript indicating that GRU.KL refers to the rate of growth 
of drug usage during the time increment from K to L. 

The nominal growth rate of addiction is determined by the number of addicts as 

NGRA.K = {■—) AD. K 7.2.3 

where 

NGRA Nominal Growth Rate of Addiction (^— 

AD Addicts (men) 

AOD Constant determining the growth rate. 

The number of drug users influences the growth rate of addiction as an availability 
multiplier of the form illustrated in Figure 7.4 where ADUM is the availability of 
drug users multiplier. 


ADUM 



Figure 7.4 Availability of Drug Users Multiplier 
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The drug education level acts to deter the flow rate and is included as a 
multiplier which decreases with increasing education effort. The form of the 
function is illustrated in Figure 7.5 where AEDM is the addiction education 
multiplier. 

AEDM 

1 .0 

2 4 6 8 10 tU 

Figure 7.5 Effect of Education Upon Addiction Growth. Rate 

The total flow from Drug Users to Addicts is the growth rate of the addiction 
level and is given as 

GRA.KL = CAEDM.K) (ADUM.K) (NGRA.K) . 7.2.4 

The population of the drug usage level is then given by 

DU. K = DU.J + (DT) (GRU. JK - GRA.JK) 7.2.5 

where GRU is the growth rate of drug usage and GRA is the growth rate of 
addiction, i.e., the flow rate from drug usage-. DT is delta time, the discrete 
time increment. 

The removal rate of addicts depends on the number of police, the effectiveness 
of police action, and the number of addicts. If it is assumed that each policeman 
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arrests a certain number of addicts per month, the nominal removal rate is given 
as 

NRRPE.K = (GAIN. K) (PE. K) 7.2.6 

where 

NRRPE Nominal Removal Rate due to Police Effort C ~ £ 2_ ) 

GAIN The effectiveness of police 

PE Police Effort (men). 

The variable "GAIN" in (7.2.6) is not a constant because addicts are increasingly 
careful as police effort increases and, as a result, police effectiveness in 
making arrests decreases. The nonlinear form of the GAIN multiplier is shown 
in Figure 7.6. 

GAIN 



Figure 7.6 Police Effectiveness 

The removal rate of addicts is also influenced by the availability of addicts to 
arrest. This effect is included as a multiplier which decreases with decreasing 
numbers of addicts, reflecting the difficulty in finding the addicts. The form 


127 



of this relationship is illustrated in Figure 7.7 where AAM is the availability 
of addicts multiplier. 


AAM 



Figure 7.7 Availability of Addicts Multiplier 


The total removal rate is then given as 

RRPE.KL = (AAM . K) (GAIN . K) (NRRPE . K) 

where RRPE is the removal rate due to police effort. The number of addicts 
then the integration of the inflow and outflow rates, i.e., 

AD.K = AD. J + (DT)( GRA.JK - RRPE. JK) . 

Police effort and the drug education program are considered to be first 
order responses to community spending. In DYNAMO this is represented as 

PE.K = PE. J + (DT) Cp^p) (CSPE . JK - PE.J) 

ED.K = ED. J + (DT)(^ig-)(CSED.JK - ED.J) 

where PE represents the Police Effort (men), DAP the Delay in Adjusting the 


7.2.7 
is 

7.2.8 

7.2.9 

7.2.10 
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Police (months), CSPE the Community Spending on Police Effort (men), CSPE the 

Community Spending on Police Effort (men), ED the Education program (men), DAP 
the Delay in Adjusting the Education program (months), and CSED the Community 
Spending on Education (men) . 

This completes the development of the system equations, however a simplifica- 
tion is now considered. The three states "Potential Users", "Drug Users", and 
"Heroin Addicts" are included in the model equations. These three states modeled 
the changing population for three divisions of the youth population. However, 
in many communities, especially those in which heroin addiction is becoming a 
problem, the time dynamics of the first two variables have been completed. That 
is, the percentage of the youth population which falls into the extremely broad 
category "Drug Users" is relatively fixed or slowly time varying, the major 
growth phase being essentially complete. For these reasons, only the variable 
"Heroin Addicts" is included as a dynamic variable. This assumption yields the 
following equations describing the system: 

i) Addicts: AD.K = AD.J + (DT) (GRA.JK - RRPE.JK) 7.2.11 

ii) Police: PE.K + PE.J + (DT) (~) (CSPE. JK - PE.J) 7.2.12 

iii) Education: ED. K = ED.J + (DT) (g^g-) (CSED. JK - ED.J). 7.2.13 


The growth rate of addiction, GRA, is given as 

GRA.KL = (AEDM.K) (Jj^)AD.K 7.2.14 

where AEDM is the effect of drug education and AOD is the nominal growth rate 
factor of addiction. The removal rate of addicts due to police effort is given 

RRPE.KL = (AAM.K) (GAIN. K) PE.K 7.2.15 
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where AAM is the availability of addicts multiplier and GAIN is the effectiveness 
of police effort. 

These discrete representations may easily be transformed into the form of 
continuous differential equations as 


f 1 C^i , x 3 ) - f^Cx^x,) = f(x 1 ,x 9 ,x 7 ) 


2 " 1’2 


!*■ 2 ’ y 


X 2 _( -DAP^ X 2 + ^DAP^ U 1 


7.2.16 


X 3 _( -DAE^ X 3 + ( -DAE' )U 2 


where x^ represents addicts, x^ police effort, x^ drug education program, u^ 
community spending on police effort, u^ community spending on drug education, 

(DAP) delay in adjusting police effort, and (DAE) the delay in adjusting the 
education program, f and f 2 represent respectively the growth rate of addiction 
and the removal rate of addicts. This system belongs to the broad class of 
nonlinear systems described as 

x = Ax + Bu + ip (x) . 7.2.17 

The results obtained in Chapters 2 and 3 regarding the optimal regulation of 
(7.2.17) will now be applied to the drug usage model. 

7.3. Optimal Regulation of the Nonlinear System 

The cost functional for the optimization problem is designed to regulate 

the number of addicts yet maintain public expenditures at a reasonable level. 

Consider the cost functional to be of the form 
T 

0 


/ 


[qxj (t) 


+ (CP)u x (t) + 


(CEju^tjjdt 


7.3.1 
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where represents addicts, community spending for police, and u 2 community 
spending for drug education. Appropriate choices for the cost parameters q, CP, 
and CE must be made to obtain "acceptable" levels of x(t) and u(t). A choice 
that is often quite reasonable [B7] is given as 


1 2 
— = maximum allowable (x,) 
q 1 J 

1 2 

Ggp) = maximum allowable (u^) 

1 2 

(^) = maximum allowable (u 2 ) 


7.3.2 


Using (2.3.8) and (2.3.9), the necessary conditions of optimality for the optimiza- 
tion problem consisting of the system (7. 2. 1,2, 3), the cost functional (7.3.1), 
and the initial condition x(0) = Xq reduce to the TPBVP 


y = Sy + if>(y) 7.3.3 

My (0) + Ny (1) = c 7.3.4 


where y is the 2n composite vector 


y = 



0 0 0 0 0 


0 


a 

(DAP) 


0 


0 - — - — 2 

(CP) (DAP) 



0 


0 


0 - a 0 

(DAP) 

0 0 0 


0 

0 


0 0 
0 0 


0 

0 


a 

(DAP) 

0 


0 

0 

a 

(CE) (DAE) 2 
0 

0 

a 

(DAE) 


7.3.5 
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a£(x 1 , x 2 , x 3 ) 


'l'(y) = 


0 

0 

-a(3f/3x 1 ) (x 1> x 2> x 3 )p 1 
-a(3f/3x 2 ) (x 1 ,x 2 ,x 3 )p 1 
-a(3f/3x 3 ) (x 1 ,x 2 ,x 3 )p 1 . 



'I 0 


0 O' 


x n] 

M = 


N = 


c = 

U 


o' 

o 


.0 I. 


o . 


7.3.6 


and where a is the change of time scale variable. 

Values must now be assigned to the various system parameters. In a sense 
these parameters depend on the community and environment being discussed. We 
assume that the community of interest is neither an extremely wealthy suburb 
nor the extremely poor section of an inner city. We assume the community has a 
population of 50,000. The youth population of such a community roughly comprises 
30% of the population [S3] . Since we are primarily interested in regulating the 
early phases of heroin usage, we assume that initially the community has a low 
level of heroin addiction, say one per thousand of the youth population. Communi- 
ties generally have a police force composed of approximately one policeman per 
thousand of population, [S3] . We assume that initially the police force has 
no effort directed specifically at heroin. The community is also assumed to 
initially have no drug education program. A reasonable value for police effective- 
ness is one conviction per month per policeman but decreasing in a nonlinear 
manner as police effort increases due to increasing caution among addicts. The 
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The delay in adjusting the police effort is essentially a training delay and is 
assumed to be six months. The delay in adjusting the drug education program is 
assumed to be one year. The boundary compatible set J = {V,M,N} must now be chosen 
for the integral representation. . The. .bound ary matrices are chosen directly from 
(7.3.7). The V matrix is chosen in the form 


ac 


ad 


ae 


0 


V = 


-aq 

0 


(DAP) 


(DAE) 


-ac 


-ad 


-ae 


(CP) (DAP) ‘ 
0 

0 


(CE) (DAE)‘ 
0 


7.3.8 


(DAP) 

0 


0 


(DAE) 


where c,d, and e may be chosen to model the nonlinearity f. The characteristic 
roots of this matrix are real, distinct, and readily evaluated, thus easing the 
determination of the P(t) matrix for convergence analysis. Numerical cases are 
now considered as examples. 

Example 7.3.9. 

In this example we consider the rather short time interval of one year. 
Specific values are selected for the cost parameters q, CP, CE, and the con- 
traction mapping method is applied to the TPBVP arising from the optimal 
regulator problem. If it is desired to prevent addiction from growing greatly 
from its initial value, q may be selected as 0.04. This represents the maximum 
desired number of addicts as 5 in (7.3.2). If the police can allocate a maximum 
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nonlinear effectiveness curve is illustrated in Figure 7.8. 

gain' 



Figure 7.8 Police Effectiveness 


The effectiveness of the drug education program is assumed to reduce the addiction 
growth rate by a maximum of 50% for a highly effective education program. The 
effectiveness is modeled as a function of the number of people involved in the 
drug education program. This is illustrated in Figure 7.9. 


AEDM 



Figure 7.9 Effect of Education Upon Addiction Growth Rate 
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of two men to the control of addiction, CP may be chosen as 0.25. Similarly, if 
the school committee believes that ten teachers are sufficient for the drug 
education program, CE may be chosen as CE = 0.01. 

The contraction mapping algorithm is begun with (t)c; y^, the center of 
S(y 0 ,r) , is chosen as the third member of the CM sequence, and r is set as 
r = 0.2. The center of SCy^.r) is shown in Figure 7.10. 



Figure 7.10 The Function y^Ct) for T = 12 Months. 
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Now determing the convergence parameters n and a , the vectors z^ and z 
composed of the elements 


2n 2n 



z i = S “P { ]C l s ij (t) + - v ± j Ct) | 

S j-1 

are evaluated as 


'0.12 ‘ 


‘0.13 ' 

0 


0 

0 

z = 

0 

0.01 


0.02 

0.11 


0.14 

.0.02. 


0.03 


Using the distinct characteristic roots, (4.7.17) is evaluated for P(t) yielding 
conservative values for n and a as 

n = sup { P (t ) z~ } = 0.14 

t 7.3.12 

a = sup {P(t)z} =0.16 . 

t 
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We then have 


r^~ = 0.17 < r = 0.2 . 

1-a 

Hence the conditions of the theorem are satisfied and the theoretical application 
of contraction mappings is successful and convergence of the CM sequence is 
indicated. The practical application of the CM algorithm reduced the convergence 

_3 

norm to 10 in ten iterations. The time histories for the state variables 
addicts, police effort, and drug education are illustrated in Figure 7.11. 



Figure 7.11 Addicts, Police, and Education for T = 12 Months 

We shall delay discussing the implications of these results until the next, example 
is presented. 

Example 7.3.13. 

In this example we consider longer term behavior and let the time interval of 
interest be four years. If it is desired to prevent addiction from growing over 
20 in the four year period, q may be selected as 0.0025. If the police can 
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allocate only one man to the control of addiction, CP may be chosen as 1.0. 
Similarly, if the school committee believes that ten teachers are sufficient 
for the drug education program, CE may be selected as 0.01. 

The contraction mapping algorithm is begun with (t)c; y Q , the center of 
S(y 0 ,r), is chosen as the fifth member of the CM sequence; and r is set as 
r = 0.2. The center of Sfy^.r) is illustrated in Figure 7.12. 



Figure 7.12 Addicts, Police, and Education for T = 12 Months 
Using (7.3.10) and S(y Q ,r), the vectors z Q and z are evaluated as 


' 0.55 ' 


' 0.59 ■ 

0 


0 

0 

z = 

0 

0.06 


0.08 

0.52 


0.54 

0.09 _ 


. 0.11 _ 
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Using a = 48, (4.7.17) is evaluated for P(t) yielding conservative values for 
n and a as 


n = sup {P(t)z } = 0.62 

t 7. 3. IS 

a = sup ,{P(t)z) = 0.73 ■ - 

t 

We have a < 1.0, however 

= 2.5 > r = 0.2 7.3.16 

1-a 

so the theoretical application of contraction mappings does not guarantee con- 
vergence. However, these results are only guidelines for the practical application 
of the CM algorithm. In fact, the CM algorithm reduced the convergence norm to 
10 in twelve iterations. The time histories for the state variables addicts, 
police effort, and drug education are illustrated in Figure 7.13, 

PE 



Figure 7.13 Addicts, Police, and Education for T = 48 Months 
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Discussion of Results 


Although the stated purpose of this chapter is illustrative in nature, 
perhaps one or two broad qualitative implicatons may be drawn from the results 
of the two sample cases. First, the need for prompt action is clearly indicated. 
With addiction growing at an exponential rate, any delay in dealing with the 
problem is critical. In both examples, the police effort with a short reaction 
time is used to begin removing the addiction core as quickly as possible. In 
the first example, it is seen that the controller responds to the short term 
situation with basically only a police effort. This is primarily due to the 
fact that the controller does not have the time to establish a viable drug 
education program. The second example is a longer term situation and the control 
response is seen to be reasonably balanced, i.e., the optimal regulator responds 
with both police effort and an education program. Again the police effort is 
the first to be utilized, but the education program is brought into play as 
quickly as possible and tends to deter long term growth. The drawing of 
quantitative conclusions from these examples would be of dubious value. However, 
the chapter illustrates that system modeling and optimal control theory may be 
jointly utilized to obtain information and insight into policy formulation for 
complex systems. Moreover, the chapter demonstrates that contraction mappings 
is a useful concept and tool for both the theoretical and practical investigation 
of nonlinear system control . 
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CHAPTER 8 


SUMMARY, CONTRIBUTIONS AND RECOMMENDATIONS 

8.1. Summary 

In the broadest sense, the objective of this dissertation was to study the 
theoretical and applied aspects of contraction mappings for the solution of 
nonlinear control problems. This objective was achieved by considering the 
theoretical and practical application of contraction mappings to the particular 
issues of optimal regulation and controllability of nonlinear dynamical systems. . 

It was shown in the study that application of the Pontryagin principle to 
the optimal regulator problem yielded necessary conditions for optimality in 
the form of a two point boundary value problem. Optimal system regulation was 
considered for both unconstrained and bounded controls and results were derived 
for the optimal regulation of linear dynamical systems and several classes of 
nonlinear systems. By an appropriate selection of boundary conditions, it was 
shown that the issue of controllability for dynamical systems may also be reduced 
to the study of two point boundary value problems. ..... 1 

The representation of two point boundary value problems by an integral 
equation was then introduced and made it possible to consider the solution of 
two point ' boundary value problems as the solution of corresponding operator 
equations. The joint application of the integral representation and the implicit 
function theorem provided new insight into the controllability of nonlinear 
systems. The methods of contraction mappings and modified contraction mappings 
were then presented for the solution of operator equations. Convergence theorems 
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were presented for both methods, and translated convergence theorems were derived 
for those operators arising from the optimal regulation of nonlinear systems. 

A detailed investigation of the calculation of the theoretical convergence criteria 
was conducted. Upper bounds were presented for the Lipschitz norm and derivative 
norm, and various techniques for evaluating these bounds were introduced. In 
particular, the use of simply structured matrices and similarity transformations 
were considered. The use of partitioned matrices in these developments provide 
considerable insight into the generic structure of the Green's matrices contained 
within the integral representation. 

Several numerical examples were presented to illustrate the theoretical and 
practical application of contraction mappings to the regulation and control of 
nonlinear systems. In particular, an example involving the regulation of 
Van der Pol's equation was used to illustrate the calculation of the convergence 
parameters and to demonstrate the manner in which the modified contraction mappings 
method may be used to extend the range of applicability of contraction mappings. 

An example considering the null controllability of the pitch motion of a satellite 
with bounded control thrust was then presented. This example illustrated the 
application of contraction mappings to an operator which did not satisfy differ- 
entiability conditions. The Lipschitz norm rather than the derivative norm was 
then used for the theoretical convergence analysis and to prove null controllabil- 
ity from the initial point. The final example involved the development of a 
dynamic model attempting to represent the causal, feedback structure of community 
drug usage. Optimal regulator theory and contraction mappings were then used 
to gain insight into how a community might best respond to a rapidly growing 
heroin addiction problem. The various examples demonstrate that contraction 
mappings is a useful tool for both the theoretical and practical investigation 
of nonlinear system control. 
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8.2. Contributions 


The author considers the following items to constitute the original contribu- 
tions of this dissertation. 

1. The determination of Green's functions in explicit form using simply 
structured matrices and similarity transformations. 

2. The development of insight into the generic structure of broad classes 
of Green's functions by the use of partitioned matrices. 

3. The development of a controllability theory for nonlinear dynamical 
systems based on an integral representation of TPBVP's, the implicit function 
theorem, and contraction mappings. 

4. The theoretical and practical application of contraction mappings to a 
nonlinear control problem with bounded input control and the subsequent use of the 
Lipschitz norm to prove convergence for the nondifferentiable operator equation. 

5. The theoretical and practical application of contraction mappings to the 
optimal regulation of a dynamic model of a socio-economic system. 

In addition, convergence theorems are presented for operators arising from 
the optimal regulation of several classes of nonlinear systems. However, these 
results are translations of the general theorems presented in Falb [FI] and in 
that sense are not completely original. 

8.3. Recommendations 

In this section some areas of possible future research will be briefly 
outlined. As indicated in the summary, the main thrust of this dissertation has 
been directed toward the application of contraction mappings. However, Falb and 
de Jong [FI] have succinctly revealed the close relationship which exists between 
contraction mappings, modified contraction mappings, and Newton's method. The 
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first area for possible additional research lies in exploiting this relationship 
and applying Newton's method to those operators arising from the optimal regula- 
tion of nonlinear systems. Investigation of the convergence criteria for Newton's 
method should yield additional insight into the theory of state regulation for 
nonlinear systems. The second area of research lies in the extension of the 
controllability results of Chapter 5. These results for the controllability of 
nonlinear systems are essentially local in nature, i.e., they consider controll- 
ability near the origin. However with additional analysis using the integral 
representation, it should be possible to identify classes of problems for which 
global results may be proved. The third and final area of recommended research 
involves an in-depth analysis into the relationship between the drug system 
model and the results of optimization. In particular, the data base for the 
model, parameter identification, and a sensitivity analysis deserve significant 
attention. In this manner, critical issues of the problem may be identified for 
additional social investigation and data collection. 


144 



APPENDIX A 


The contraction mappings program consists of a main program and several 
subroutines. A brief description of the function of each part is now presented. 

MAIN essentially directs the program and performs no actual computation. 

MAIN first calls the subroutine STTRM which calculates the fundamental matrices 
<t^(t,0) and <J>^(0,s). To accomplish this task, STTRM calls AFCT and VELEMS, 
and the integration to calculate is performed by DIFEQ. The resultant 

fundamental matrices are stored by OUTP and OUTT. MAIN next calls CALC, the 
major subroutine of the algorithm. CALC computes the Green's functions and 
directs the solution of the successive members of the CM sequence. VCAL and 
VELEMS are used to calculate V(t), and SBFN calculates (c - g[y(0)) - h(y(0)) + 
My(0) + Ny (1) and F(y)}. FINT then calls DQSF to integrate the expression 

t 1 

/ Gj[t ,s){F(y n (s) ,s)-V(s)y n (s) }ds . /ej, (t,s){F(y n (s),s)-V(s)y n (s)}ds . 
0 t 

CONV is then called to test for convergence. If the test for convergence is 
successful, the program returns to MAIN and ends. If the test for convergence 
fails, the algorithm remains in CALC and calculates the successive solutions 
until either convergence is attained or a stop condition is reached. All 
computations are done in double precision arithmetic. To use the contraction 
mappings program, the user must modify only two subroutines, VELEMS and SBFN. 

In VELEMS, the user specifies the choice of the V(t) matrix. In SBFN, 
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the user specifies the differential equation y = F(y,t) and the boundary 
condition g(y(0)) + h(y(l)) = c. The program contains many comment statements 
to ease application. 
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r» o o 


CONTRACTION MAPPING ALGORITHM 
M A T N 

DOUBLE. P P F C I S 1 0 M PH I , PH I S , PH I US » DEL T , FN » D , YS » 3 I NT , QQ I NT , V , C 
UOUBL E PPECrSIDM XN,XM 
DOUBLE PRECISION FMINT 
DOUBLE PRECISION UNITY 
DOUBLE PRECISION UN 

COMMON PH I ( 9,0,? I ),PHIS(9, Q ,21) » DEL I , FN ( 9 » 2 1 ) , 0 ( 9 ) 

COMMON Y S ( 9 » 2 1 » 1 5 ) , Q I N T ( 9 » 2 1 ) » OQ l NT ( 9 » 2 1 ) » V ( 9 , p , 2 1 ) 

COMMON G I 9 ) , XN I 9,9) ,XM( 9 ,9) , 1 1 , I I I 
COMMON kK,LL,NDM,NINT, ITFR 
DIMENSION UNI TY ( 15 , IS ) 

C N D I M IS THE DIMENSION OF THE PROBLEM VECTOR 
RE AO (E, 2 ) NDTM 

2 FORMAT ( IB) 

NDM=NDTM 

N SO = ND I M* NO I M 

C THE INCREMENT OF SOLUTION IS NOW READ IN. 

READ IS, 3) OFLT 

3 c 0 R M A T ( D 1 0 . ? ) 

C THF NUMBER OF INCREMENTS IS NOW CALCULATED. 

F N I N T = 1 . 0 DO /DEL T+ 1 . 100 
NINT=IDI NT{ FNINT) 

C THE SUB ROUT IMF STTRM WILL NOW BE CALLED TO CALCULATE THE STATE 
C TRANSITION MATRIX (IF THF SPECIFIED LINEAR SYSTEM AND ITS ADJOIN 
CALI STTRMINDFM) 

C THF MATRIX UNITY IS FORMED TO CHECK THE ACCURACY IN CALCULATING, 
C PHI AND PHIS. 

DO 66 3 J= 1 , MOM ■ - 

DO 663 1=1, NDM 
UN l T Y ( I , J )=0.000 
DO 663 KM, NDM 

66 3 UNITYU, J) = UNITY! I , J ) +PH I ( I , K , 2 1 I * D H I S ( K , J , 2 11 
DO 665 1=1, NDM 

WRTTE(6,666) (UNITY(T,J), J= 1 , NDM ) 



664 FORMAT ( * *,5X,015.9) 

665 CONTINUE 
UN-0 .ODO 

DO 3 33 1=1 , NO l M 
UN= UNITY! I« I )HJN 
333 CONTINUE 

IF (UN .GT. 1.5 00* NO I M ) GO TO 606 

C NOW THF MAJOR SUBROUTINE CALC IS CALLED TO CALCULATE AND STUFF T H r: 
C NFW SOLUTION. 

CALL CALC(NDIM) 

C A STOP CONDITION IS CHECKED. 

IF ( ITER . EQ. 15) Gil TO 6 06 
ITER l = ITCR+ 1 
DO 19 K=1,ITER1 
W R I T F ( 6 , 1 .6 ) K 

14 FORMAT ( 'OS 6X,4HK = ,13) 

DO IB J = 1 , N 0 1 M 
WRITE (6,2 0) J 

30 c 0 P M A T ( » • , 10X,4HJ = ,13) 

WP I T E ( 6 , 1 7 ) ( YS ( J , NOS t K ) , NDS=1,NINT) 

1 7 FORMAT! • • , 15X ,D1 5. 8) 

19 CONTINUE 
19 CONTINUE 
606 STOP 
END 
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r> o 


SUBROUT I Ni- SBFN(NDIM) 


THIS SUBROUTINE CALCULATE S T H P VALUES OF FN = F ( Y ) -V *Y FUR VALUES OF MOT 
DOUBLE PR FT is I • J M PH T , PH I S , PH I OS , DEL T , FN, D , YS , 0 I NT , QQ I NT , V ,C 
DOUBLE PRECISION XN,XM 

DOUBLE PRECISION Y V , F , Z , Y I * Y T , G , H f T T , TT T 

DOUBLE PRECISION X 1 , X? T X 3 , X/- , X5 , X6 , DFOX T , G A I M , GMAX.DGDP 

OOUBL C PR PC T SION GT &U , AO 0 , &L DM, P I , AAM, A VB t DR DA , PCS , F F 
DOUBLE PRECISION C .A P , D AF , CP f CF , NI , 02 , Q3 
COMMON PH I ( 9,9 , 21 ) f PHIS ( <3,9,2 I ) ,OELT, FN( 9,21 ) ,D (9 ) 

COMMON YS ( 9, 2L , IS) ,QTNT ( 9 , 21 ) ,00 1 NT (9,21) , V ( 9 , 9 , 2 1 ) 

C OMMON Cl 9 ) , XN ( 9 , 9 ) , XM ( 9 , 0 ) , 1 1 , I I I 
COMMON K K , L L , NO M , N I N T , I T F P. 

D I MF; MS T on Y V U. 5 ) , F ( 1 5 ) , Z H 5 ) , Y T ( 1 5 ) , Y T ( l 5 ) , G ( 1 5 ) , H ( I 5 ) 

D l M= MS I ON TT( 15 ) , TTT (1 5 ) 

dimension dfdx r (9,4 ) ,pcs( a » ,ff( qi 

C THE FOLLOWING VARIABLES A R c USED TO CALCULATE THE NONLINEAR 
C FHKC.TNG FUNCTION F(I). 
n l =3. 141 59300 
AVB=40 . DDO 
DT=1?.000 
GMAX= 1.000 
GTAlJs; 25500 
DA p-6, 000 
DAE-12.0D0 
C E = . 0100 
C°-.?500 
01 =.0400 
02=0.000 
03=0.0 D 0 

no 6 0 0 NO T = 1 , N I NT 

C THF NONLINEAR FOUATION IS \ FUNCTION OF THE STAFF AT THE CURRENT 
C TIME. A VECTOR OF THr STAFF AT THF NOT IS C.REAIFO AMD IS IJSFD TO 
C CALCULATE F AT NOT. 

00.59 9 T=l,NDIM 
599 YV ( I ) =YS ( I .NOT , ITER ) 
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x l = y v ( l ) 

X 2 = Y V ( 2 ) 

X3=YV ( 3 ) 

X4=YV (4 ) 

X5 = YV ( 5 ) 

X6= YV ( 6 ) 

00 102 1-1*3 
00 102 J = 2, 3 
102 0F0XT ( I » J )=0.0D0 

G A I N= DC OS ( 0 • 3600* X2 ) 

OGO P = -0 • 36D0*0S I N ( 0 . 3600 * X? ) 

A0D=5 0. 000 

TF(X3 .GT. 10.0 DO) GO TO 110 
AEOM=. 7500+.25DO*OCbS(PI *X3/10.000) 

OFOXT ( 3 , 1 )-=-( • 2 50 0*P I/10.000)*DSIN( P l *X3 / 1 0. 000 ) *X 1 / AOD 
GO TO 111 
UO AFOM = , 5 000 

OFDXT ( 3 » 1 ) = 0 • ODD 

111 C ONT I NUF 

IF (XI .GT. 60.000) GO TO 112 
AAM-.5D0+.5Q0*DSIM(PI*{X1-AV8/2.0D0)/4VB> 
DR0A=GAIN*X2*.5DO*(PI/AVB)*DCOS(PI*(Xl-AVB/2.000)/AVB) 
GO TO 113 

112 AAM=1.000 
DRD A=0. 000 

113 CONTINUE 

OFDXT ( 1 ,1 ) = AFDM/ AOD-DROA 
OFOXT (2 ,1 )=-AAM*GAIN-AAM*X?*OGOP 
DO 109 K = 1 1 3 
10« PCS(K)=YV(K+3) 

00 l OR 1=1,3 
FF ( I )=0.000 
DO 1 C 8 K = 1 , 3 

108 FF(I)=FF( I ) +OFOXT ( I ,K)*PCS(K) 

GRA=AEDM*X1/A0D 

RRPE=AAM*GAIN*X2 
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F< l)=DT*M GRA-RRPF) 

F(2 )=DT*(-.l ,0D0*XR/DAP-1 . 000*X5/t CAP*0AP*CP) 1 
P( $ ) =PT*( - X3 /DA E-X6 / { CAE*DAE*CE ) ) 

F(4)=QT*(-QI*X1-FF( 1)1 
F ( 5 )=DT*( -Q?*X2*X5/DAP-FF (?) ) 

F(6 )=DT*( -Q3*X3+-X6/DAF-EFt T) ) 

C THF PRODUCT V ( NOT 1 *Y ( NOT ) TS NOW 03T Af NED AMO STORED AS A 
C FUNCTION OF TTMF. 

00 598 1= ltNOIM 
U I )=0.0D0 
on 598 K= ltNOIM 
59 8 U I ) = 7 ( I )4-V( I v K«Nf)T)*YV( K) 

C FN IS MOW OBTAINED AND STORED AS A FUNCTION OF TIME 
DO 59 7 Tl = UNDIM 
FP7 FN( TL,NOT)=F( IL 1-Z( II 1 
C THIS PROCEDURE IS REPEATED FOR INCREASING NOT 
500 CONTINUE 

C THIS SUBROUTINE ALSO CALCULATES THE EXPRESSION, 

C 0=C-Gm-H(Y)+XM*Y(0)+XN*Y(l). THE INITIAL AND TERMINAL STATE 
C VECTORS ARE GENERATED BELOW. 

00 601 I = l » NOT M 
Y I ( I ) = Y S ( I, 1 , ITER) 

601 Y T ( I ) =YS< I , N lNT t ITER ) 

G(l) = YI ( l ) 

G ( 2 ) - Y I ( R ) 

g m * y i m 

G ( A ) =0. ODD 
G ( 5 ) = C . ODO 
G ( 6 ) =0. 000 
H ( 1 ) =0. 000 
H (2) =0.000 
M ( 3 ) =0.000 
H ( 4 ) = YT ( 4 ) 

H(5)=YT(5» 

H ( 6 ) = YT ( 6 J 

C THE PPOOUCTS M * Y ( 0 ) AND N*Y(1) ARE NOW OBTAINED AND THF RESULT n 



C IS FORMED. 

^0 602 1=1,NPIM 
TT( I) =0.000 
00 602 K=i,NDIM 

60 2 TT(T )=TT( II+XM( I,K)*YT(K) 
nn 603 1=1* NOT M 

TTT ( I >=0.000 
on 603 K = 1 , N D 1 

603 TTT<n = TTTtI) «-XN( I ,K) *YT(K) 

00 604 M=l»NOIM 

604 D ( M ) = C ( M ) -G ( M > -H ( M ) +TT ( M ) +TTT ( M ) 
RETURN 

FNO 



SUBROUTINE VELEMS(X,A) 

C THIS SUBROUTINE CALCULATES THE LINEAR SYSTEM MATRIX IN VECTOR FORM 
DOUBLE PRECISION X,A 
DOUBLE PRECISION TF 

DOUBLE PRECISION A 1 , A2 , A 3 , DT , C AP., C A E , Q 1 , 02 , Q3, C P , C F 
DIMENSION A (22 5) 

A 1 =.0200 
A 2 = - . I 5 D 0 
A3=-.0O5D0 
DT = 12 .000 
QAP=6. ODO 
OA f = l 2 . ODO' 

C E = • 0 IDO 
CP=.25D0 
91 = . 0400 
0?=C. ODO 
03 - 0 . one 

A ( I ) = A].*QT 
A ( ? ) = C . 000 
A( 31*0.000 
A ( A ) = - 0 1 w D T 
A ( 5 ) = C . Of) 0 
M 0 ) = 0. ODO 
A ( 7 l=A2*0T 
A ( H ) =-:)T/OAP 
.A (0 ) - 0 . 0 D 0 
A ( 1 01 -0. ODO 
A ( .1 1 )=-Q2*DT 
A ( 12 ) *0.000 
A ( ]. 3 ) = A 3 * 0 T 
A ( 14 )=C.000 
A ( ] 5 ) = - 0 T / D A F 
A ( 16 ) =0.0 DO 
A (17) *0.000 
A ( 1 fil =-03*0 T 
A(] 9 1=0. ODD 
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A ( 20 ) =0.000 
A ( 2 1 ) =0. 000 
A(22)=-Al*nT 
A(?3 )=•- A2*DT 
A ( 24 ) =- A3 *1.)T 
A(25>=0.000 

A126)=-DT/(DAP*CAP*CP) 

A { 27 ) =0.000 

A ( 28 )=0.0D0 

A(?9)=0T/DAP 

A ( 30 = 0.000 

A(71)=0.000 

A (32) =0.000 

A ( 3 3 >=-0T /( 0 A E *DA E *C F ) 

A ( 34 ) =0. 000 

A ( 35 ) =0.000 

A ( 36 ) = 0T / DA E 

PFUJRN 

PND 
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SUP ROUTINE VC AL 

C THIS SUBROUTINE CALCULATES and STORES THE V MATRIX IN TIME. 
DOUBLE PRECISION X,A 

DOUBLE PRECISION PH I , PH l S , °H I OS , DEL T , FN , 0 , Y S , 0 I NT , QQ I N T , V , C 
DOUBLE PRECISION XN,XM 

COMMON PH I ( 9,9, 21 > ,PHIS( 9,9,21 ) ,0ELT,P.N(9,21 ) , 0(9 ) 

COMMON YS( 9, 21 , 1 S) ,Q INTI Q,21 ) ,OClNT (9,21 ) ,V( 9,9,21) 

COMMON C (9 ) , XN( 9, 9) , XM( 9 ,9 ) , 1 1 , I I T 
COMMON KK,LL ,NDM,NINT, ITEP 
DIMENSION A ( 2 2 5 ) 

DO 100 J = 1 , 2 1 
x=( J-D^OELT 
CALL VELEMS(X,A) 

DO 555 rO=l,NDM 
DO 554 J0= 1 , NOM 
K Q= { I Q- 1 ) *NOM+JQ 
559 VI JO, 10, J ) = A(KQ) 

559 CGNTINUF 
100 CONTINUE 
RETURN 
ENO 
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SUBROUTINE AFCT ( X,SM,DERV) 

C THIS SUBROUT INF IS USED to CALCULATE V( . ) IN THE INTEGRATION 
C FOR PHI . 

DOUBLE PRECISION PH l , P HI S, PH I OS, DEL T , F N , D , YS , 3 I NT , 001 NT , V , C 

DOUBLE PRECISION XN, XM 

DOUBLE PRECISION X , SM , DE RV , A , TM° 

COMMON PH I (9,9, ?1 ) ,PH IS ( 9,9,21) ,DELT,FN( 9,21) ,0(9) 

COMMON YS (9,21, 15) ,QINT( 9,21 ) ,QQTNT( 9,21 ) ,V(9, 9, 21) 

COMMON C ( 9 ) ,XN(9»9),XM(9,9) ,11, I II 
COMMON KK, LL ,NDM,N TNT, ITER 
DIMENSION DERVU5) ,A(??5 ),SM(15) 

C A I L VELEMSI X, A ) 

DO 5 5 5 I Q = 1 , N DM 
TMP=0. ODD 
DO 554- JQ = 1 , NDM 
K0= ( JO-l) *NDM+TQ 
554 TMP=TMP+A (KQ) *SM< JO) 

5*5 D E P V ( 1 0 ) - T M P 
RETURN 
END 
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SUBROUTINE ATFC ( X , SM ,DERV ) 

C THIS SUBROUTINE IS USED TO CALCULATE -V • ( . > IN THF INTEGRATION 
C FOR PHIS. 

DOUBLE PRECISION PHI , PHI S, PH I OS » DEL T , FN , D ♦ YS ,Q l NT , QQ I NT , V ,C 

DOUBLE PRECISION XN,XM 

nnUBl.F PRECISION X , SM , OF P V , A , T MP 

COMMON PH I (9,9, 21 » ,PHIS( 9,9,21 J *DEL T,FN(9,2l» t 0(9) 

COMMON VS(9,?1 ,15> ,'0INT( 9,2] ),QQINT( 9, 21) ,V(9, 9, 21) 

COMMON C ( 9 ) » XN ( 9 . 0) ,XM( Q ,9) , 1 1 , I T T 
COMMON K.K,LL,NDM,NINT, JTFR 
0 I ME NS I ON OF RV U 5 ) , A < .? 25 ) , SM ( 1 5 > 

CALL VELEMSIX, A I 
DO 555 10=1, WDM 

tmp=o. ono 
no 554 J0= 1 , NOM 
K 0= ( IQ-1 ) *NDM+JQ 

554 TMP=TMP+A(KQI*SM( JQ) 

555 DEE V ( IQ )= -TMP 
RFTUEN 

END 
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SUBROUTINE D I FEQ ( N, PMODE , T ,DT ,CTR ,VAR,RHS ) 
c this subroutine is used TO INTEGRATE FOR PHI AMO PHIS. 

C THF TECHNIQUE IS A FOURTH ORDER RUNGE-KUTTA AS MODIFIED BY GILL. 
DOUBLE PRECISION V AR ( 6 ) , RHS ( 2 ) , QL AM ( 50 \ , CC.C1 , C CC2 ,CCC B 
DOUBLE PRECISION UGHL Y f R 00T2 , MNUS , PLUS 
DOUBLE PRECISION T,DT 

30 FORMAT (43HIMP ROPER COUNTER SETTING IN THE DIFFQ SURRO) 
INTEGER CTR, PMODE 
I F ( P MODE ) 9 0,1,2 
l DO 4 J= 1 , N 

4 QL AM ( J ) -0 . 

R 00 T 2 = 1 • 4 1 4 2 1 3 5 6 2 3 7 3 0 ° 5 0 0 
MNUS=1. D0-l.no/R00T2 
PLUS -1 . DO+1 .DO/ROOT? 

PM0DE=1 
CTR = 0 

? IF(CTR) 99,3,5 
3 CCC I = . EDO 
CCC 2= 1 . DO 
CCC3 =DT* . 5D0 
T = T+CCC 3 
GO TO 20 

5 I F ( C TR-2 I 6,7,3 

6 CCC 1= MNUS 
14 CC.C2 = CCC1 

CCC 3=CCC 1 T 
GO TO 20 

7 CCC 1 =PLUS 
T=T+OT* , 5D0 
GO TO 14 

P CCC 1=. 166666666666666700 
CCC 2= . 333333333333333300 
CCC3=DT*.5D0 
C TR =- 1 

20 CTR=CTR+I 
CCC 1=CCC1*DT 



DO 22 J=l *N 

•JGHL Y = CCC 1*PHS( J> -CGR 2*QLAM( J ) 

OL AM { Jt -QI.AM( J ) +IJGHL Y +UGHLY+UGHL Y-CCC3*R HS ( J ) 
2? VAR( J)=VAR( JJ-KJGHI.Y 
RETURN 

°9 WRITF<6,T0) 

RETURN 

END 
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SUBROUTINE STTRM(NDIM) 

r. this subroutine computes the state transition matrix or the linear 

C SYSTEM AND ITS ADJOINT AND STORES THEM AS FUNCTIONS GF TIME. 

OOUBL E PR EC I S I ON PH I , PH I S , PH I OS , DEL T , FN, 0 , YS , Q I NT f QQ I NT , V , C 

DOUBLE PRECISION XN,XM 

DOUBLE PRECISION Y , DE: R Y , TF , T , DT 

C OM MON PH l ( 9 , 9 , 2 1 ) , PH I S ( 9 , 9 , 2 1) , DEL T , FN ( 9 , 2 U , D < 9 ) 

COMMON YS(9 t 21 , 151 ,QINT( <5,21 ) ,OQINT(9 t 21 ) »V(9,9, 21) 

COMMON f. (9) » XN ( 9, 9) » XM ( 9 ,9 ) , II. I IT 
• COMMON KK,LL ,NDM,NTNT,TTfP 
DIMENSION Y(] 5) ,DERY( 15) 

INTEGER CTR,PMODE 
READ (5,1) OT 

1 FORMAT (010. 2 ) 

WRITE (6,2) DT 

2 FORMAT ( '0 5X , 19HINTEGRATT0N STEP = ,015.9) 

00 7 11=1, NOIM 

00 3 J = I , NO I M 

1 F ( ( I I- J ) . EQ. 0) Y( J) =1.000 
I ^ ( < I T - J ) . N E . 0) Y ( J ) =0 . 000 

3 CONTINUE 
T=0 . ODO 
TF = 1 .000 
PMOQE^o 
CT R =0 

KK = 0 

CALL OU TP ( T , Y , NO I M ) 

A CONTINUE 

CALL AFCT(T, Y,OERY) 

CALL DIFEQ(NOIM,PMODE,T,DT,CTR»Y,UERY) 

I F ( CTR .EQ. 0) GO TO 5 
GO TO 4 

5 CALL OUTP (T,Y,NOIM) 

TF (T .GE. TF) GO TO 6 
GO TO 4 

6 continue ■ ■; 
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7 CONT TNUF 

DO 9 NN=1 » ND I M 

WRITF(6»8) IPHIINN,NM f ?l), NM=1,NDIM) 

8 FORMAT ( 'O' ,5X, D15.5} 

9 CONT TNUF 

DO l A I I 1 = 1, NOT M 
00 10 J = 1 » ND I M 

IF((HI-.n .FQ. 0 ) Y C J ) = 1 • 000 
T F ( ( I I I - J ) • NE . 0) Y ( J ) = 0 . ODO 
10 CONTINUE 
T= 0.000 
T c = 1 .000 
PMOOF-O 
CTR = 0 
LL = 0 

CALL OUTT (T,Y f NDlM) 

It CONT I NUF 

CALL A T F C ( T t Y » 0 E P Y ) 

CALL DIFEO( N D I M » PMOOE » T t DT »CTR , Y » DE R Y ) 
IFICTP ,PQ. 0) GO TO 12 
CO TO 11 

12 CALL OUTTI T,Y,MDIM) 

I F ( T -,GE . TFI GO TO 13 
GO TO 1 1 

13 CONTINUE 
19 CONTINUE 

DO 16 NN= 1 1 NDI M 

WRT TF ( 6 1 1 8 ) (PHIS (NN,NMt2l >, NM=l»NCIM) 

15 FORMAT! 'O' , 5X.015.8) 

16 CONTINUE 
RETURN 
END 
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SUBROUTINE OUT P (T,Y,NDIM) 

C THIS SUBROUTINE STORES THE MATRIX PHII.,0! AT THF apppopriatf 
C INCREMENTS OF TIME. 

DOUBLE PR EC I SION PH I , PHI S , PH I OS , DEL T , F N, 0, Y S , 0 I NT , 001 NT , V , 

OOUBLF PR C C IS TON XN,XM 

DOUBLE PRFCISION T,Y 

DOUBLE PRECISION DELT , QQ, TEST, DABS 

COMMON P4I(9,9,21),PHIS(9,9,21),DELT,FNI9,21),D(9) 

COMMON YS < 9 , 21 , 15 ) ,Q INT ( 9 ,21) ,QQINT ( 9,? 1 ) ,V (9, 9, 21 ) 

COMMON C(9),XN<Q,9t,XM(9,9),Il,lIl 
COMMON KK,LL, NOM,NTNT, ITER 
DIMENSION Y ( 9 ) 

OELT = .0500 
QQ = F LOAT ( KK ) 

TEST-DABS ( QQ*DELT-T) 

I F ( TES T .GT. .000100) GO TO 100 
WRT TE (6,1 01 ) T 
101 FOR MAT (• ' , 4HT = ,015.8) 

KK = KK +1 

DO 99 J = 1 , N D I M 
PHI(J,TT,KK)=Y(J) 

Q9 CONTINUE 
100 CONTINUE 
RETURN 
FNO 



SUBROUTINE OiJTT ( T » Y , ND I M > 

C THIS SUBROUTINE STORES THE MATRIX PHTS(.,0) AT APPROPRIATE 
.C INCREMENTS OF TIME. 

OOURL.E PRECISION PH I , P H I S, PH I OS , DEL T , FN , D, YS , 0 1 NT, QQ I NT , V ,C 

OOUBLE PRECISION XN,XM 

DOUBLE PRECISION T,Y 

DOUBLE PRECISION DELT ,RR , TEST , DAB S 

COMMON PHI (9, 9,21 1,PHIS( 9,9,2 U, DEL T,FN< 9,211, 0<Q) 

COMMON YSI9.21, 151 fOTNT(9,?l|,QQiNT(9,2l| ,v(9,9,21l 
COMMON C(9) v XN(9 v 9)«XM(9«9) f II, ITT 
COMMON K K , L L , NDM , M I N T , ITER 
DIMENSION Y ( 9 ) 

DELT=.05DO 
B R = E L O.A TILL) 

TEST = DABS (RR*DELT-T) 

IE (TEST . GT. .0001 DO). GO TO 100 
WRITFl6,i01) T 
101 FORMAT ( ' ' , 4HT = ,015. R) 

LL=LL+l 

00 99 J = 1 , ND I M 
PHIS! Ill , J,LL»=Y< J ) 

99 CONTINUE 
LOO CONTINUE 
RFTURN 
END 
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SUB ROUT INF CALCCNOIM) 

C THTS IS THF MAJOR SUBROUTINE TN THE PROGRAM. HERE THE INTEGRAL 
C E OUA T TONS ARE SOLVED FOR THE ITERATED SOLUTIONS AND THE TFST FOR 
C CONVERGENCE IS MADE. 

DO’HJL £ 

DOUBLE 
DOUR l E 
DOUBL E 
DOUBLE 
DOUBLE 
DOUBL E 
DOUBL E 
DOUBL F 
COMMON 
COMMON 
C OMMON 
C 0 m m 0 N 

DIMENSION TEMPI ( 15), TEMP? ( 15>,TEMPR( 1 5 ) , L ( 1 5 ) , M ( 1 5 ) 

DIMENSION T ( 1 5 , lf> , 21 ) , TM ( l 5 , l 5 , .?1 ) , TN< 1 5 » 1 5 , 2 1 > 

DIMENSION XC ( 15 , 15 ) ,VS I ( 725) , SI I N( I 5,15 ) ,SI ( 15 , 15) 

DIMENSION TSI ( ? ? 5 ) ,TSIIN(2?5) 

DIMENSION 0 1 V ( 2 1 ) ,Z(2l) 

FQIJ I V A|. ENCE ( S I ( 1 , 1 ) , TSI ( 1 ) ) 

EQUI V AL ENCF ( S I I N ( l , 1 ) , T S I IN { 1 ) ) 

C THE RELAXATION FACTOR IS READ IN. NORMALLY TT IS ONE. 

PFAH(5,555) RE 
555 c ORMAT(Dl 0.2 ) 

C THE BOUNDARY CONDITION MATRICES IN THE BOUNDARY COMPATIBLE 
C SET J = ( V , M , N ) ARE NOW READ IN. 

DO 2 I - 1 , NO I M 

RE AD ( 5,1 ) ( XM ( I , J ) , J = l. ,NDIM) 

1 FORMAT (010. 2) 

? CONTINUE 

DO A 1 = 1 , NO I M 

RF AD ( 5 , B ) ( X N ( I , J ) , J = 1 , NO I M ) 

? FORM AT (DIO. 2) 


RP EC IS I DM PHI, PHI S, PH l OS, DEL T , EM , D, YS , iil NT ,001 NT ,V ,C 
PRECISION X N , X M 

PR EC I SI ON S I , XC ,T ,TM,TN, TEMPI ,TFMP2 , TEMPO 

PRECISION EPS 

PRECISION DET 

PRECISION VS I, SUN 

PRECISION TSI, TSI IN 

PRECISION RF 

PRECISION COST, CIV, 7 

PH I ( P,9, ,?1 ) , PHIS ( 9,0 ,21 ) ,DELT ,FN(R,2l ) , D(9 ) 

YS ( o, ?1 , 15) ,Q INT( o , 21 ) ,OCINT( 9,21 ) ,V( 9, 9,21 ) 

C (9) ,XN(9,9) ,XM(c t q» t n, II I 
KK,LL ,NDM,NINT f ITER 
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4 CONTINUE 

C EPS. THE CONVERGENCE MEASURE IS NOW READ IN. 

RF.A0(5,P) EPS 
? P 0 P M A T (DIO..?) 

C THE ROUND AP Y CONDITION VECTOR 0 IS NOW READ IN. 

P E A D ( 5 , 10 ) (C(I)r I = 1 » ND I M ) 

10 FORMAT (01 0.2) 

C THE CONSTANT ISM IS NOW READ IN. IF ISM IS ONE, THE PROGRAM 

C COMPUTES THE INITIAL BOUNDARY COMPATIBLE GUESS. 

f. IF ISM is NOT ONE, THE INITIAL SOLUTION- IS NOW PEA 0 IN. 

READ (5, 666) ISM 
66 6 POP MATH 10) 

C NOW FORMING THE PRODUCT OF N*PHI(1,0) 

DO 7 J = 1 , N 0 T M 
DO 7 1=1. ND I M 
XC( I, J 1=0.000 
DO 7 K = 1 , ND I M 

7 XC ( I » J ) =XC ( r . J ) ♦- X N ( I ,K)* PHI ( K, J ,NINT ) 

C THE MATRIX SUM ( M+N*PH I ( 1 , 0) ) ) IS NOW FORMED. 

DO P J = 1 , NO I M 
DO 5 1 = 1, NO I M 

a s i ( r , ji = xm( i , j h-xc< i , j ) 

WRT TE (6,12 ) 

12 FORMAT ( *0 ', 2X f 2HS I ) 

DO 15 1=1. ND I M 

DO I** J = 1 . M D I M 
WRITE (6,13) S I ( I,J) 

13 FORMAT ( ■ • ,10X,01 5.8) 

14 CONTINUE 

15 CONTINUE 
mode = 7 

CALL APR. A Y ( MODE , ND I M , NO I v , 1 5 , 1 5 , V S I , T S I ) 

CALL M l n v ( VS I , N D I M » D E T , L , M ) 

M0DF=1 

CALL ARRAY(MODE,NDIM,NOIM f 15,15, VSI ,TSI IN) 

WP.I TE (6, 1 12 ) 
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11? F ORMAT I ' 0 ' , 2X «4HS 1 1 N ) 

DO It. 5 I = 1 » N D I M 
00 L 1 A J=1,NDIM 
WRI TE (6,1 13) SIIN(I,J) 

M3 FORMAT! 'O' , 10X, 01 5. 8) 

114 CONTINUE 

11 5 CONTINUE 

00 4CO N0T=1,NINT 
00 397 J=l,NOIM 
DO 397 1=1, N DIM 
T ( I , J , NOT ) =0. ODO 
DO 39 7 K = 1 * NO I M 

39 7 T ( I , J , NOT ) = T ( I , J , NOT ) + PH I ( I , K , NO.T ) * S 1 1 N I K , J ) 
00 398 J = l » N D I M 
00 398 I =1 , NO l M 

TM ( I ♦ J , NOT ) = 0 . ono 

00 398 K=1,NDIM 

7 98 TM ( I , J , N T T ) = T 3 ( T , J , NO T ) + T ( l , K , NOT ) * X.V { K , J » 

00 399 J=l, NO I M 
00 399 1=1, NO I M 
TNI I ♦ J,NDT) =0.000 

00 399 K=l,NOIM 

7 09 TNI I , J , NO T ) = T N I I , J , NOT J + T I I , K ,NOT )*XC( K, J ) 

400 CONTINUE 
I T ER = 0 

401 I TER = I TFR + 1 

1 f= I ITER . EO. 13 ) GO TO 9 09 
TF< ITER . GT . 1 ) GO TO 73 2 

T F I ISM .EQ. 1 ! GO to 669 
0.0 668 I=l,NOIM 

RF AO I 5 » 6671 (YS(I,J,1), J = 1,NINT) 

667 FORMATIOIO. 2) 

668 CONTINUE 
GO TO 670 

669 DO 814 NDS= 1 ,N I NT 
00 813 T = 1 , N D I M 
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YS( I ,N0S, 1 ) =0.000 
DO 816 K= 1 1 MO I M 

816 YS< I ,NDS, 1 ) =YS< I ,NDS , l)+T{ I ,K,MDS)*C(K) 

PI 5 CONTINUE 
814 CONTINUE 
670 CONTINUE 

00 818 1 = 1, NDI M 

WRITE (6» 817 ) ( YS( I ,N, 1 ) , N=1,NINT) 

81 7 FORMAT ( • • , 20X , Dl 5 . 8 1 

818 CONTINUE 

THE SUBROUTINE VCAL IS MOW CALLED TO CALCULATE AND STORE THE LINFAP 
SYSTEM MATRIX AS A FUNCTION OF TIME 
CALL VCAL 

SUBROUTINE S BEN WILL MOW BE CALLED TO CALCULATE FN=F(Y|-V*Y 
782 CALL SBFN(NDIM) 

SUBROUTINE PINT INTEGRATES P H I ( 0 ♦ S ) *F N ( S ) FROM ZERO TO T AND 
STORES THE INTEGRAL AS A FUNCTION OF T, WHERE T VARIES FROM ZERO 
C TO ONE. THESE VALUES ARE USFO TCi CALCULATE THE INTEGRA! FROM T TO ONE 
CALL FINT(NDIM) 

C THE NEXT SEQUENCE OF INSTRUCTIONS SOLVES FOR THE NEXT ITERATED 
C SOLUTION. FIRST THE PRODUCT T ( T J *D WILL BE CALCULATED. 

DO 304 NDS = I » N I NT 
DO 300 I = l , N D I M 
T EM P 1 ( I J - 0 , 0 DO 
DO 300 K= 1 , ND I M 

30C TFMPK I ) = TEMPI! I )+T( I , K , NDS > *D ( K ) 

C NEXT, THE PPOOUCT Qt TM(NDS) AND THE INTEGRAL OF FN FROM ZERO TO 
C NOS IS FORMED. 

DO 301 T= l » NO I M 
TFMP2 (I )=0.0D0 
DO 301 K= 1 , NO I M 

301 TFMP2(!> = TEMP2( I) *TM< I t K,NOS)*QINT(K,NDS) : 

C NEXT, THE PRODUCT OF TNINDSJ AND THE INTEGRAL OF FN FROM NDS TO ONE 
C IS FORMED. 

DO 302 1=1 , NDI M • 

TFMP.3 { I ) =0. ODO 



\ 


DO 302 K= 1 » NO T M 

30? T E M P 3 ( I ) = TFMP3( I) +TN { I » K » NOS ) *QQ I NT ( K , NDS ) 

C THF- THREE TEMPS ARE SUMMED TO GIVE THE VALUE OF THE NFW SDLUTTUM 
C AT TIME NOS. 

DO 303 J J = 1 » iMD I M 

303 YS( JJ ,NDS, I TFR 1 1) = ( 1 . ODO-R F > * YS ( J J , NDS , ITFR ) 

C ♦RF’rs (TEMPI ( JJ ) +TFMP2( JJ ) - T E M P 3 ( J J ) ) 

r. THTS PR Of. ED UR F TS REPEATED FOR INCREASING NDS 

304 CONTINUE 

I T E R 1 - T T F R + 1 
WRITE (6, A 04 ) 

404 F0OMAT ( »0 ' , 5X,2HYS» 

DO 407 I = l t ND I M 

DO ^06 N- 1 , N f NT 

WRITE (6i405I YS ( I , N ♦ I TFR 1 ) 

405 FORMAT!' • , 1 OX , D1 5 . 8 ) 

406 CONTINUE 

407 CONTINUE 

C CONVERGENCE OF THF ITERATION IS NOW TESTED. 

CALL CONV (N DIM, MM, EPS) 

IF (MM . FQ. 1 ) GO T 0 401 
909 RETURN 
END • 


168 



SUBROUTINE FINT(NDfM) 

C SUBROUTINE PINT INTEGRATES PHI ( 0» S) *FN( S) FROM ZERO TO T AND 
C STORES THE INTEGRAL AS A FUNCTION OF T, WHERE T VARIFS FROM 7 FRO 
C TO ONE. THESE VALUES ARF USED TO CALCULATE THE INTEGRAL FROM T TO ONE. 
nPUPI F PRECISION f>HI t PHIS,PMlOS,DELT t FN,D,YS,OIMT t OOlNT f V t C 
DCURI E PR F C I SION XN.XM 
DOUBLE PRECISION OC,7,CTV 

COMMON PHI (9,9, 21 ) »PHIS( 9,9,2] > , D FI T , FN ( 9 , ? 1 ) .0(9) 

COMMON YS( 9 ,21 , 15 ) ,0 INTI 9 ,21 ) ,00 I NT (9,2 1 ) »V(9, 9, 21 ) 

COMMON C(9),XN(9,9),XM(9,9),II,I II • 

COMMON KK.I. L.NDM.NINT, IT C R 
DIMENSION 0 0 ( 1 5 . ? 1 ) » Z ( 21 ) » 0 1 V ( 2 1 ) 

C CALCULATE AND STORE THE VECTOR PH I ( 0 . NOS ) *FN ( NOS ) AS A FUNCTION OF NDS 
DO 97 N D S = 1 , N J N T 
90 95 I = 1 , N D I M 
OQ ( I , NDS) =0.000 
DO 94 K = l , N 0 I M 

94 QQ ( I , NDS ) =00 ( I , NDS ) ■‘•PH f S ( I , K , NO S ) *F N ( K , NOS ) 

95 CONTINUE 

97 CONTINUE 

C THE TIME HISTORY OF FACH COMPONENT IS PUT IS VECTOR FORM AND 
C INTEGRATED BY OQSF. 

00 100 !<J=1,M0TM 
00 Qfl LJ=1,NINT 

0 I V ( L J ) =DQ ( K J , L J ) 

98 CONTINUE 

CALL DOSF ( OF LT , 01 V t Z » N IN T ) 

C THE INTEGRALS ARE STORED IN CT NT AND 00 1 NT . 

DO 99 NN= 1 , N INT 

0 1 N T ( K J . N N ) = Z ( M N ) 

99 CONTINUE 
100 CONTINUE 

DO 202 M-l.NDIM 
DO 20 1 MM= 1 , N IN T 

201 QOINT (m,mm) = 01 NT (M, NT NT )-OINT( M.MM) 

202 CONTINUE 
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SUBROUTINE CONV ( NOT M f MM, FP S ) 

C TUTS SUBROUTINE TESTS FOR CONVERGENCE OF THE ITERATION. 

OOlJRI. E PRECISION PH I , PHI S , PH I OS , OEL T , F N, p, YS , 0 I NT , 00 INT , V , C 

DOUBLE PRECISION XN,XM 

DOUBLE PRECISION OY , B IGC. , B TG , EPS , CON 

OOUBLF PRECISION DABS 

DOUBLE PRECISION C0ST,QIV,7 

COMMON PH I (9,9, 21 ), PH IS (9,9,21 ) t DEL T , FN ( 9 , 2 1 ) ,0(9) 

COMMON YS<9 , 21 , 15) ,QINT( 9,21 ) ,QOINT( 9,21 > ,V( 9,9,21 ) 

C. DM MON C ( 9 > , XN ( 9 , 9 ) , XM ( 9 , 9 ) , 1 1 , I I I 
COMMON KK,L L,NDM,NINT, ITER 
DIMENSION 0 Y ( 2 I ) ,9IGC( 1 5 ) 

DIMENSION 0IVI21) ,7121) 
on 700 I = l , N 0 1 M 
no 693 ND S= 1 , N I NT 

6 98 0Y(NDS)=DA8S(YS( I , NOS , ITFP+1)-YS( I, NOS, ITER) ) 

C THE LARGEST ABSOLUTE DIFFERENCE IN THIS COMPONENT WILL NOW V 
C THE LARGEST ABSOLUTE DIFFERENCE IN THIS COMPONENT WILL NOW BE FOUND 
B I G = D Y ( 1 ) 

DO 699 N= 2, N IN T 
T F ( D Y ( M ) .LT. BIG) GO TO 699 
B I G = 0 Y ( M ) 

699 COMTINUF 

700 BIGCI I ) = B I G 
CON = BIGr. I 1 ) 

00 7 C 1 1=2, NO I M 

1 F ( B I GC ( I. ) .LT. CON) GO TO 701 
C 0N = B I GC ( L ) 

701 CONTINUE 

WRITF (6,755) CON 

7^5 FORMAT! *0', 15X, 30HN0RM OF FUNCTION DIFFERENCE = ,015.8) 

IF (CON .LT. EPS) GO TO ° 9 <* 

.MM = 1 

GO TO 998 

999 

998 RETURN 
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